From 6fce6d0245e1b37feb8a032a114ef7655c0150cb Mon Sep 17 00:00:00 2001
From: Lukas Beiske <lukas.beiske@tu-dortmund.de>
Date: Mon, 30 Sep 2024 15:43:59 +0200
Subject: [PATCH 1/2] Make things compatible with numpy 2.0; replace deprecated
 logging.warn()

---
 pyirf/benchmarks/energy_bias_resolution.py  | 11 +++--
 pyirf/binning.py                            | 36 ++++++++------
 pyirf/compat.py                             | 10 ++++
 pyirf/cuts.py                               | 26 +++++++----
 pyirf/interpolation/component_estimators.py | 15 ++++--
 pyirf/sensitivity.py                        | 52 ++++++++++++---------
 pyirf/spectral.py                           |  9 ++--
 pyirf/statistics.py                         |  5 +-
 pyirf/utils.py                              |  3 +-
 setup.py                                    |  3 +-
 10 files changed, 107 insertions(+), 63 deletions(-)
 create mode 100644 pyirf/compat.py

diff --git a/pyirf/benchmarks/energy_bias_resolution.py b/pyirf/benchmarks/energy_bias_resolution.py
index 675c9af18..4479549c3 100644
--- a/pyirf/benchmarks/energy_bias_resolution.py
+++ b/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)
diff --git a/pyirf/binning.py b/pyirf/binning.py
index 1aa3ecd02..0eb4413c4 100644
--- a/pyirf/binning.py
+++ b/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,
     )
 
diff --git a/pyirf/compat.py b/pyirf/compat.py
new file mode 100644
index 000000000..a52cd60a5
--- /dev/null
+++ b/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
diff --git a/pyirf/cuts.py b/pyirf/cuts.py
index 4f83268b1..ffda0e906 100644
--- a/pyirf/cuts.py
+++ b/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)
diff --git a/pyirf/interpolation/component_estimators.py b/pyirf/interpolation/component_estimators.py
index cd14beade..06e123f37 100644
--- a/pyirf/interpolation/component_estimators.py
+++ b/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,
         )
diff --git a/pyirf/sensitivity.py b/pyirf/sensitivity.py
index a8f7b114d..c8873c828 100644
--- a/pyirf/sensitivity.py
+++ b/pyirf/sensitivity.py
@@ -1,6 +1,7 @@
 """
 Functions to calculate sensitivity
 """
+
 import numpy as np
 from scipy.optimize import brentq
 import logging
@@ -8,6 +9,7 @@
 from astropy.table import QTable
 import astropy.units as u
 
+from .compat import COPY_IF_NEEDED
 from .statistics import li_ma_significance
 from .utils import check_histograms, cone_solid_angle
 from .binning import create_histogram_table, bin_center
@@ -36,7 +38,7 @@ def _relative_sensitivity(
         return np.nan
 
     if n_on < 0 or n_off < 0:
-        raise ValueError(f'n_on and n_off must be positive, got {n_on}, {n_off}')
+        raise ValueError(f"n_on and n_off must be positive, got {n_on}, {n_off}")
 
     n_background = n_off * alpha
     n_signal = n_on - n_background
@@ -65,7 +67,7 @@ def equation(relative_flux):
         relative_flux = brentq(equation, lower_bound, upper_bound)
 
     except (RuntimeError, ValueError) as e:
-        log.warn(
+        log.warning(
             "Could not calculate relative significance for"
             f" n_signal={n_signal:.1f}, n_off={n_off:.1f}, returning nan {e}"
         )
@@ -87,8 +89,7 @@ def equation(relative_flux):
 
 
 _relative_sensitivity_vectorized = np.vectorize(
-    _relative_sensitivity,
-    excluded=['significance_function']
+    _relative_sensitivity, excluded=["significance_function"]
 )
 
 
@@ -230,8 +231,10 @@ def calculate_sensitivity(
 
     # convert any quantities to arrays,
     # since quantitites don't work with vectorize
-    n_on = u.Quantity(n_on, copy=False).to_value(u.one)
-    n_off = u.Quantity(background_hist["n_weighted"], copy=False).to_value(u.one)
+    n_on = u.Quantity(n_on, copy=COPY_IF_NEEDED).to_value(u.one)
+    n_off = u.Quantity(background_hist["n_weighted"], copy=COPY_IF_NEEDED).to_value(
+        u.one
+    )
 
     # calculate sensitivity in each bin
     rel_sens = relative_sensitivity(
@@ -257,7 +260,9 @@ def calculate_sensitivity(
     s["n_background_weighted"] = background_hist["n_weighted"]
 
     # copy also "n_proton" / "n_electron_weighted" etc. if available
-    for k in filter(lambda c: c.startswith('n_') and c != 'n_weighted', background_hist.colnames):
+    for k in filter(
+        lambda c: c.startswith("n_") and c != "n_weighted", background_hist.colnames
+    ):
         s[k] = background_hist[k]
 
     s["significance"] = significance_function(
@@ -271,10 +276,14 @@ def calculate_sensitivity(
 
 
 def estimate_background(
-    events, reco_energy_bins, theta_cuts, alpha,
-    fov_offset_min, fov_offset_max, 
+    events,
+    reco_energy_bins,
+    theta_cuts,
+    alpha,
+    fov_offset_min,
+    fov_offset_max,
 ):
-    '''
+    """
     Estimate the number of background events for a point-like sensitivity.
 
     Due to limited statistics, it is often not possible to just apply the same
@@ -308,33 +317,30 @@ def estimate_background(
         Minimum distance from the fov center for background events to be taken into account
     fov_offset_max: astropy.units.Quantity[angle]
         Maximum distance from the fov center for background events to be taken into account
-    '''
-    in_ring = (
-        (events['reco_source_fov_offset'] >= fov_offset_min)
-        & (events['reco_source_fov_offset'] < fov_offset_max)
+    """
+    in_ring = (events["reco_source_fov_offset"] >= fov_offset_min) & (
+        events["reco_source_fov_offset"] < fov_offset_max
     )
 
     bg = create_histogram_table(
         events[in_ring],
         reco_energy_bins,
-        key='reco_energy',
+        key="reco_energy",
     )
 
     # scale number of background events according to the on region size
     # background radius and alpha
     center = bin_center(reco_energy_bins)
     # interpolate the theta cut to the bins used here
-    theta_cuts_bg_bins = np.interp(
-        center,
-        theta_cuts['center'],
-        theta_cuts['cut']
-    )
+    theta_cuts_bg_bins = np.interp(center, theta_cuts["center"], theta_cuts["cut"])
 
-    solid_angle_on = cone_solid_angle(theta_cuts_bg_bins) 
-    solid_angle_ring = cone_solid_angle(fov_offset_max) - cone_solid_angle(fov_offset_min)
+    solid_angle_on = cone_solid_angle(theta_cuts_bg_bins)
+    solid_angle_ring = cone_solid_angle(fov_offset_max) - cone_solid_angle(
+        fov_offset_min
+    )
     size_ratio = (solid_angle_on / solid_angle_ring).to_value(u.one)
 
-    for key in filter(lambda col: col.startswith('n'), bg.colnames):
+    for key in filter(lambda col: col.startswith("n"), bg.colnames):
         # *= not possible due to upcast from int to float
         bg[key] = bg[key] * size_ratio / alpha
 
diff --git a/pyirf/spectral.py b/pyirf/spectral.py
index 35bad8647..182657bff 100644
--- a/pyirf/spectral.py
+++ b/pyirf/spectral.py
@@ -1,6 +1,7 @@
 """
 Functions and classes for calculating spectral weights
 """
+
 import sys
 from importlib.resources import as_file, files
 
@@ -9,6 +10,7 @@
 from astropy.table import QTable
 from scipy.interpolate import interp1d
 
+from .compat import COPY_IF_NEEDED
 from .utils import cone_solid_angle
 
 #: Unit of a point source flux
@@ -120,7 +122,9 @@ def from_simulation(cls, simulated_event_info, obstime, e_ref=1 * u.TeV):
         viewcone_max = simulated_event_info.viewcone_max
 
         if (viewcone_max - viewcone_min).value > 0:
-            solid_angle = cone_solid_angle(viewcone_max) - cone_solid_angle(viewcone_min)
+            solid_angle = cone_solid_angle(viewcone_max) - cone_solid_angle(
+                viewcone_min
+            )
             unit = DIFFUSE_FLUX_UNIT
         else:
             solid_angle = 1
@@ -308,7 +312,6 @@ def __init__(
         self.interp = interp1d(x, y, bounds_error=False, fill_value="extrapolate")
 
     def __call__(self, energy):
-
         x = (energy / self.reference_energy).to_value(u.one)
 
         if self.log_energy:
@@ -319,7 +322,7 @@ def __call__(self, energy):
         if self.log_flux:
             y = 10**y
 
-        return u.Quantity(y, self.flux_unit, copy=False)
+        return u.Quantity(y, self.flux_unit, copy=COPY_IF_NEEDED)
 
     @classmethod
     def from_table(
diff --git a/pyirf/statistics.py b/pyirf/statistics.py
index b4f5b10de..988bfd70b 100644
--- a/pyirf/statistics.py
+++ b/pyirf/statistics.py
@@ -1,5 +1,6 @@
 import numpy as np
 
+from .compat import COPY_IF_NEEDED
 from .utils import is_scalar
 
 __all__ = ["li_ma_significance"]
@@ -34,8 +35,8 @@ def li_ma_significance(n_on, n_off, alpha=0.2):
 
     # Cast everything into float64 to avoid numeric instabilties
     # when multiplying very small and very big numbers to get t1 and t2
-    n_on = np.array(n_on, copy=False, ndmin=1, dtype=np.float64)
-    n_off = np.array(n_off, copy=False, ndmin=1, dtype=np.float64)
+    n_on = np.array(n_on, copy=COPY_IF_NEEDED, ndmin=1, dtype=np.float64)
+    n_off = np.array(n_off, copy=COPY_IF_NEEDED, ndmin=1, dtype=np.float64)
     alpha = np.float64(alpha)
 
     with np.errstate(divide="ignore", invalid="ignore"):
diff --git a/pyirf/utils.py b/pyirf/utils.py
index d5ed730e5..28e1da5eb 100644
--- a/pyirf/utils.py
+++ b/pyirf/utils.py
@@ -2,6 +2,7 @@
 import astropy.units as u
 from astropy.coordinates import angular_separation
 
+from .compat import COPY_IF_NEEDED
 from .exceptions import MissingColumns, WrongColumnUnit
 
 
@@ -27,7 +28,7 @@ def is_scalar(val):
     result: bool
         True is if input object is a scalar, False otherwise.
     """
-    result = np.array(val, copy=False).shape == tuple()
+    result = np.array(val, copy=COPY_IF_NEEDED).shape == tuple()
     return result
 
 
diff --git a/setup.py b/setup.py
index 561acc80d..e4ad60003 100644
--- a/setup.py
+++ b/setup.py
@@ -41,12 +41,13 @@
 
 setup(
     use_scm_version={"write_to": os.path.join("pyirf", "_version.py")},
-    packages=find_packages(exclude=['pyirf._dev_version']),
+    packages=find_packages(exclude=["pyirf._dev_version"]),
     install_requires=[
         "astropy>=5.3,<7.0.0a0",
         "numpy>=1.21",
         "scipy",
         "tqdm",
+        "packaging",
     ],
     include_package_data=True,
     extras_require=extras_require,

From 49f2b4e131fdd5e82dbd5f55034928d512a43382 Mon Sep 17 00:00:00 2001
From: Lukas Beiske <lukas.beiske@tu-dortmund.de>
Date: Mon, 30 Sep 2024 15:53:05 +0200
Subject: [PATCH 2/2] Add changelog entry

---
 docs/changes/289.maintenance.rst | 3 +++
 1 file changed, 3 insertions(+)
 create mode 100644 docs/changes/289.maintenance.rst

diff --git a/docs/changes/289.maintenance.rst b/docs/changes/289.maintenance.rst
new file mode 100644
index 000000000..379a7e43f
--- /dev/null
+++ b/docs/changes/289.maintenance.rst
@@ -0,0 +1,3 @@
+Make pyirf compatible with numpy 2.0.
+
+Replace deprecated ``logging.warn`` with ``logging.warning``.