Merge pull request #289 from cta-observatory/numpy2_compat

Add compatibility with numpy 2.0 and replace deprecated logging.warn()

docs/changes/289.maintenance.rst

@@ -0,0 +1,3 @@
+Make pyirf compatible with numpy 2.0.
+
+Replace deprecated ``logging.warn`` with ``logging.warning``.

pyirf/benchmarks/energy_bias_resolution.py

@@ -4,6 +4,7 @@
 import astropy.units as u
 
 from ..binning import calculate_bin_indices, UNDERFLOW_INDEX, OVERFLOW_INDEX
+from ..compat import COPY_IF_NEEDED
 
 
 NORM_LOWER_SIGMA, NORM_UPPER_SIGMA = norm(0, 1).cdf([-1, 1])
@@ -83,8 +84,10 @@ def energy_bias_resolution(
     """
 
     # create a table to make use of groupby operations
-    table = QTable(events[["true_energy", "reco_energy"]], copy=False)
-    table["rel_error"] = (events["reco_energy"] / events["true_energy"]).to_value(u.one) - 1
+    table = QTable(events[["true_energy", "reco_energy"]], copy=COPY_IF_NEEDED)
+    table["rel_error"] = (events["reco_energy"] / events["true_energy"]).to_value(
+        u.one
+    ) - 1
 
     energy_key = f"{energy_type}_energy"
 
@@ -102,7 +105,6 @@ def energy_bias_resolution(
         # we return the table filled with NaNs
         return result
 
-
     # use groupby operations to calculate the percentile in each bin
     bin_index, valid = calculate_bin_indices(table[energy_key], energy_bins)
     by_bin = table.group_by(bin_index)
@@ -115,6 +117,7 @@ def energy_bias_resolution(
         result["resolution"][bin_idx] = resolution_function(group["rel_error"])
     return result
 
+
 def energy_bias_resolution_from_energy_dispersion(
     energy_dispersion,
     migration_bins,
@@ -147,7 +150,7 @@ def energy_bias_resolution_from_energy_dispersion(
             low, median, high = np.interp(
                 [NORM_LOWER_SIGMA, MEDIAN, NORM_UPPER_SIGMA],
                 cdf[energy_bin, :, fov_bin],
-                migration_bins[1:] # cdf is defined at upper bin edge
+                migration_bins[1:],  # cdf is defined at upper bin edge
             )
             bias[energy_bin, fov_bin] = median - 1
             resolution[energy_bin, fov_bin] = 0.5 * (high - low)

pyirf/binning.py

@@ -7,6 +7,8 @@
 import astropy.units as u
 from astropy.table import QTable
 
+from .compat import COPY_IF_NEEDED
+
 
 #: Index returned by `calculate_bin_indices` for underflown values
 UNDERFLOW_INDEX = np.iinfo(np.int64).min
@@ -37,12 +39,12 @@ def join_bin_lo_hi(bin_lo, bin_hi):
         The joint bins
     """
 
-    if np.allclose(bin_lo[...,1:], bin_hi[...,:-1], rtol=1.e-5):
-        last_axis=len(bin_lo.shape)-1
-        bins = np.concatenate((bin_lo, bin_hi[...,-1:]), axis=last_axis)
+    if np.allclose(bin_lo[..., 1:], bin_hi[..., :-1], rtol=1.0e-5):
+        last_axis = len(bin_lo.shape) - 1
+        bins = np.concatenate((bin_lo, bin_hi[..., -1:]), axis=last_axis)
         return bins
     else:
-        raise ValueError('Not matching bin edges')
+        raise ValueError("Not matching bin edges")
 
 
 def split_bin_lo_hi(bins):
@@ -62,10 +64,11 @@ def split_bin_lo_hi(bins):
     bin_hi: np.array or u.Quantity
         Hi bin edges array
     """
-    bin_lo=bins[...,:-1]
-    bin_hi=bins[...,1:]
+    bin_lo = bins[..., :-1]
+    bin_hi = bins[..., 1:]
     return bin_lo, bin_hi
 
+
 def add_overflow_bins(bins, positive=True):
     """
     Add under and overflow bins to a bin array.
@@ -124,7 +127,7 @@ def create_bins_per_decade(e_min, e_max, bins_per_decade=5):
     # include endpoint if reasonably close
     eps = step / 10000
     bins = 10 ** np.arange(log_lower, log_upper + eps, step)
-    return u.Quantity(bins, e_min.unit, copy=False)
+    return u.Quantity(bins, e_min.unit, copy=COPY_IF_NEEDED)
 
 
 def calculate_bin_indices(data, bins):
@@ -156,7 +159,7 @@ def calculate_bin_indices(data, bins):
 
     valid: np.ndarray[bool]
         Boolean mask indicating if a given value belongs into one of the defined bins.
-        False indicates that an entry fell into the over- or underflow bins. 
+        False indicates that an entry fell into the over- or underflow bins.
     """
 
     if hasattr(data, "unit"):
@@ -169,8 +172,8 @@ def calculate_bin_indices(data, bins):
     n_bins = len(bins) - 1
     idx = np.digitize(data, bins) - 1
 
-    underflow = (idx < 0)
-    overflow = (idx >= n_bins)
+    underflow = idx < 0
+    overflow = idx >= n_bins
     idx[underflow] = UNDERFLOW_INDEX
     idx[overflow] = OVERFLOW_INDEX
     valid = ~underflow & ~overflow
@@ -211,11 +214,11 @@ def create_histogram_table(events, bins, key="reco_energy"):
         hist["n_weighted"] = hist["n"]
 
     # create counts per particle type, only works if there is at least 1 event
-    if 'particle_type' in events.colnames and len(events) > 0:
-        by_particle = events.group_by('particle_type')
+    if "particle_type" in events.colnames and len(events) > 0:
+        by_particle = events.group_by("particle_type")
 
         for group_key, group in zip(by_particle.groups.keys, by_particle.groups):
-            particle = group_key['particle_type']
+            particle = group_key["particle_type"]
 
             hist["n_" + particle], _ = np.histogram(group[key], bins)
 
@@ -282,8 +285,11 @@ def resample_histogram1d(data, old_edges, new_edges, axis=0):
     cumsum /= norm
 
     f_integral = interp1d(
-        old_edges, cumsum, bounds_error=False,
-        fill_value=(0, 1), kind="quadratic",
+        old_edges,
+        cumsum,
+        bounds_error=False,
+        fill_value=(0, 1),
+        kind="quadratic",
         axis=axis,
     )
 

pyirf/compat.py

@@ -0,0 +1,10 @@
+import numpy as np
+from packaging.version import Version
+
+
+# in numpy 1.x, copy=False allows copying if it cannot be avoided
+# in numpy 2.0, copy=False raises an error when the copy cannot be avoided
+# copy=None is a new option in numpy 2.0 for the previous behavior of copy=False
+COPY_IF_NEEDED = None
+if Version(np.__version__) < Version("2.0.0.dev"):
+    COPY_IF_NEEDED = False

pyirf/cuts.py

@@ -4,17 +4,25 @@
 import astropy.units as u
 
 from .binning import calculate_bin_indices, bin_center
+from .compat import COPY_IF_NEEDED
 
 __all__ = [
-    'calculate_percentile_cut',
-    'evaluate_binned_cut',
-    'compare_irf_cuts',
+    "calculate_percentile_cut",
+    "evaluate_binned_cut",
+    "compare_irf_cuts",
 ]
 
 
 def calculate_percentile_cut(
-    values, bin_values, bins, fill_value, percentile=68, min_value=None, max_value=None,
-    smoothing=None, min_events=10,
+    values,
+    bin_values,
+    bins,
+    fill_value,
+    percentile=68,
+    min_value=None,
+    max_value=None,
+    smoothing=None,
+    min_events=10,
 ):
     """
     Calculate cuts as the percentile of a given quantity in bins of another
@@ -46,7 +54,7 @@ def calculate_percentile_cut(
     """
     # create a table to make use of groupby operations
     # we use a normal table here to avoid astropy/astropy#13840
-    table = Table({"values": values}, copy=False)
+    table = Table({"values": values}, copy=COPY_IF_NEEDED)
     unit = table["values"].unit
 
     # make sure units match
@@ -84,10 +92,10 @@ def calculate_percentile_cut(
             cut_table["cut"].value[bin_idx] = value
 
     if smoothing is not None:
-        cut_table['cut'].value[:] = gaussian_filter1d(
+        cut_table["cut"].value[:] = gaussian_filter1d(
             cut_table["cut"].value,
             smoothing,
-            mode='nearest',
+            mode="nearest",
         )
 
     return cut_table
@@ -126,7 +134,7 @@ def evaluate_binned_cut(values, bin_values, cut_table, op):
         passes the bin specific cut given in cut table.
     """
     if not isinstance(cut_table, QTable):
-        raise ValueError('cut_table needs to be an astropy.table.QTable')
+        raise ValueError("cut_table needs to be an astropy.table.QTable")
 
     bins = np.append(cut_table["low"], cut_table["high"][-1])
     bin_index, valid = calculate_bin_indices(bin_values, bins)

pyirf/interpolation/component_estimators.py

@@ -1,16 +1,19 @@
 """Classes to estimate (interpolate/extrapolate) actual IRF HDUs"""
+
 import warnings
 
 import astropy.units as u
 import numpy as np
 from scipy.spatial import Delaunay
 
+
 from .base_extrapolators import DiscretePDFExtrapolator, ParametrizedExtrapolator
 from .base_interpolators import (
     DiscretePDFInterpolator,
     ParametrizedInterpolator,
     PDFNormalization,
 )
+from ..compat import COPY_IF_NEEDED
 from .griddata_interpolator import GridDataInterpolator
 from .nearest_neighbor_searcher import BaseNearestNeighborSearcher
 from .quantile_interpolator import QuantileInterpolator
@@ -428,9 +431,9 @@ def __init__(
         self.min_effective_area = min_effective_area
 
         # remove zeros and log it
-        effective_area[
-            effective_area < self.min_effective_area
-        ] = self.min_effective_area
+        effective_area[effective_area < self.min_effective_area] = (
+            self.min_effective_area
+        )
         effective_area = np.log(effective_area)
 
         super().__init__(
@@ -468,7 +471,7 @@ def __call__(self, target_point):
         # remove entries manipulated by min_effective_area
         aeff_interp[aeff_interp < self.min_effective_area] = 0
 
-        return u.Quantity(aeff_interp, u.m**2, copy=False)
+        return u.Quantity(aeff_interp, u.m**2, copy=COPY_IF_NEEDED)
 
 
 class RadMaxEstimator(ParametrizedComponentEstimator):
@@ -852,5 +855,7 @@ def __call__(self, 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
+            np.swapaxes(interpolated_psf_normed, -1, self.axis),
+            u.sr**-1,
+            copy=COPY_IF_NEEDED,
         )