Add support for numpy 2.0

docs/changes/2580.feature.rst

@@ -0,0 +1 @@
+ctapipe is now compatible with numpy 2.0.

pyproject.toml

@@ -38,7 +38,7 @@ dependencies = [
     "joblib",
     "matplotlib ~=3.0",
     "numba >=0.56",
-    "numpy ~=1.16",
+    "numpy >=1.23,<3.0.0a0",
     "psutil",
     "pyyaml >=5.1",
     "requests",
@@ -48,6 +48,7 @@ dependencies = [
     "tqdm >=4.32",
     "traitlets ~=5.6",
     "zstandard",
+    "packaging",
 ]
 
 [project.optional-dependencies]

src/ctapipe/atmosphere.py

@@ -19,6 +19,8 @@
 from astropy.table import QTable, Table
 from scipy.interpolate import interp1d
 
+from .compat import COPY_IF_NEEDED
+
 __all__ = [
     "AtmosphereDensityProfile",
     "ExponentialAtmosphereDensityProfile",
@@ -355,17 +357,19 @@ def __init__(self, table: Table):
     @u.quantity_input(height=u.m)
     def __call__(self, height) -> u.Quantity:
         log_density = self._density_interp(height.to_value(u.km))
-        return u.Quantity(10**log_density, DENSITY_UNIT, copy=False)
+        return u.Quantity(10**log_density, DENSITY_UNIT, copy=COPY_IF_NEEDED)
 
     @u.quantity_input(height=u.m)
     def integral(self, height) -> u.Quantity:
         log_col_density = self._col_density_interp(height.to_value(u.km))
-        return u.Quantity(10**log_col_density, GRAMMAGE_UNIT, copy=False)
+        return u.Quantity(10**log_col_density, GRAMMAGE_UNIT, copy=COPY_IF_NEEDED)
 
     @u.quantity_input(overburden=u.g / (u.cm * u.cm))
     def height_from_overburden(self, overburden) -> u.Quantity:
         log_overburden = np.log10(overburden.to_value(GRAMMAGE_UNIT))
-        return u.Quantity(self._height_interp(log_overburden), u.km, copy=False)
+        return u.Quantity(
+            self._height_interp(log_overburden), u.km, copy=COPY_IF_NEEDED
+        )
 
     def __repr__(self):
         return (

src/ctapipe/compat.py

@@ -3,6 +3,9 @@
 """
 import sys
 
+import numpy as np
+from packaging.version import Version
+
 __all__ = [
     "StrEnum",
 ]
@@ -15,3 +18,11 @@
 
     class StrEnum(str, Enum):
         """Compatibility backfill of StrEnum for python < 3.11"""
+
+
+# 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

src/ctapipe/coordinates/camera_frame.py

@@ -15,6 +15,7 @@
     frame_transform_graph,
 )
 
+from ..compat import COPY_IF_NEEDED
 from .representation import PlanarRepresentation
 from .telescope_frame import TelescopeFrame
 
@@ -144,10 +145,14 @@ def camera_to_telescope(camera_coord, telescope_frame):
     # where theta is the angle to the optical axis and r is the distance
     # to the camera center in the focal plane
     fov_lat = u.Quantity(
-        (x_rotated / focal_length).to_value(u.dimensionless_unscaled), u.rad, copy=False
+        (x_rotated / focal_length).to_value(u.dimensionless_unscaled),
+        u.rad,
+        copy=COPY_IF_NEEDED,
     )
     fov_lon = u.Quantity(
-        (y_rotated / focal_length).to_value(u.dimensionless_unscaled), u.rad, copy=False
+        (y_rotated / focal_length).to_value(u.dimensionless_unscaled),
+        u.rad,
+        copy=COPY_IF_NEEDED,
     )
 
     representation = UnitSphericalRepresentation(lat=fov_lat, lon=fov_lon)

src/ctapipe/image/extractor.py

@@ -1709,7 +1709,7 @@ def __call__(
             peak_time = adaptive_centroid(
                 d_waveforms, peak_index, leading_edge_rel_descend_limit
             )
-            peak_time = peak_time / (self.sampling_rate_ghz[tel_id] * upsampling)
+            peak_time /= self.sampling_rate_ghz[tel_id] * upsampling
 
         if gain != 0:
             charge /= gain

src/ctapipe/instrument/subarray.py

@@ -16,6 +16,7 @@
 from astropy.utils import lazyproperty
 
 from .. import __version__ as CTAPIPE_VERSION
+from ..compat import COPY_IF_NEEDED
 from ..coordinates import CameraFrame, GroundFrame
 from .camera import CameraDescription, CameraGeometry, CameraReadout
 from .optics import FocalLengthKind, OpticsDescription
@@ -206,7 +207,7 @@ def tel_ids_to_indices(self, tel_ids):
         np.array:
             array of corresponding tel indices
         """
-        tel_ids = np.array(tel_ids, dtype=int, copy=False).ravel()
+        tel_ids = np.array(tel_ids, dtype=int, copy=COPY_IF_NEEDED).ravel()
         return self.tel_index_array[tel_ids]
 
     def tel_ids_to_mask(self, tel_ids):
@@ -710,7 +711,7 @@ def from_hdf(cls, path, focal_length_choice=FocalLengthKind.EFFECTIVE):
 
         optic_descriptions = [
             OpticsDescription(
-                name=row["optics_name"],
+                name=str(row["optics_name"]),
                 size_type=row["size_type"],
                 reflector_shape=row["reflector_shape"],
                 n_mirrors=row["n_mirrors"],
@@ -745,7 +746,7 @@ def from_hdf(cls, path, focal_length_choice=FocalLengthKind.EFFECTIVE):
 
             camera.geometry.frame = CameraFrame(focal_length=focal_length)
             telescope_descriptions[row["tel_id"]] = TelescopeDescription(
-                name=row["name"], optics=optics, camera=camera
+                name=str(row["name"]), optics=optics, camera=camera
             )
 
         positions = np.column_stack([layout[f"pos_{c}"] for c in "xyz"])
@@ -761,7 +762,7 @@ def from_hdf(cls, path, focal_length_choice=FocalLengthKind.EFFECTIVE):
             )
 
         return cls(
-            name=name,
+            name=str(name),
             tel_positions={
                 tel_id: pos for tel_id, pos in zip(layout["tel_id"], positions)
             },

src/ctapipe/io/simteleventsource.py

@@ -20,6 +20,7 @@
     TableAtmosphereDensityProfile,
 )
 from ..calib.camera.gainselection import GainChannel, GainSelector
+from ..compat import COPY_IF_NEEDED
 from ..containers import (
     ArrayEventContainer,
     CoordinateFrameType,
@@ -179,8 +180,12 @@ def build_camera(simtel_telescope, telescope, frame):
     )
     readout = CameraReadout(
         telescope.camera_name,
-        sampling_rate=u.Quantity(1 / pixel_settings["time_slice"], u.GHz),
-        reference_pulse_shape=pixel_settings["ref_shape"].astype("float64", copy=False),
+        sampling_rate=u.Quantity(
+            1.0 / pixel_settings["time_slice"].astype(np.float64), u.GHz
+        ),
+        reference_pulse_shape=pixel_settings["ref_shape"].astype(
+            "float64", copy=COPY_IF_NEEDED
+        ),
         reference_pulse_sample_width=u.Quantity(
             pixel_settings["ref_step"], u.ns, dtype="float64"
         ),
@@ -932,18 +937,18 @@ def _fill_event_pointing(tracking_position):
 
         # take pointing corrected position if available
         if np.isnan(azimuth_cor):
-            azimuth = u.Quantity(azimuth_raw, u.rad, copy=False)
+            azimuth = u.Quantity(azimuth_raw, u.rad, copy=COPY_IF_NEEDED)
         else:
-            azimuth = u.Quantity(azimuth_cor, u.rad, copy=False)
+            azimuth = u.Quantity(azimuth_cor, u.rad, copy=COPY_IF_NEEDED)
 
         # take pointing corrected position if available
         if np.isnan(altitude_cor):
             altitude = u.Quantity(
-                _clip_altitude_if_close(altitude_raw), u.rad, copy=False
+                _clip_altitude_if_close(altitude_raw), u.rad, copy=COPY_IF_NEEDED
             )
         else:
             altitude = u.Quantity(
-                _clip_altitude_if_close(altitude_cor), u.rad, copy=False
+                _clip_altitude_if_close(altitude_cor), u.rad, copy=COPY_IF_NEEDED
             )
 
         return TelescopePointingContainer(azimuth=azimuth, altitude=altitude)

src/ctapipe/io/tableio.py

@@ -9,6 +9,8 @@
 from astropy.time import Time
 from astropy.units import Quantity
 
+from ctapipe.compat import COPY_IF_NEEDED
+
 from ..core import Component
 from ..instrument import SubarrayDescription
 
@@ -316,7 +318,7 @@ def __call__(self, value: Time):
         return getattr(getattr(value, self.scale), self.format)
 
     def inverse(self, value):
-        return Time(value, scale=self.scale, format=self.format, copy=False)
+        return Time(value, scale=self.scale, format=self.format, copy=COPY_IF_NEEDED)
 
     def get_meta(self, colname):
         return {
@@ -336,7 +338,7 @@ def __call__(self, value):
         return value.to_value(self.unit)
 
     def inverse(self, value):
-        return Quantity(value, self.unit, copy=False)
+        return Quantity(value, self.unit, copy=COPY_IF_NEEDED)
 
     def get_meta(self, colname):
         return {
@@ -399,8 +401,8 @@ def __init__(self, scale, offset, source_dtype, target_dtype):
             self.maxval = iinfo.max - 1
 
     def __call__(self, value):
-        is_scalar = np.array(value, copy=False).shape == ()
-        value = np.atleast_1d(value).astype(self.source_dtype, copy=False)
+        is_scalar = np.array(value, copy=COPY_IF_NEEDED).shape == ()
+        value = np.atleast_1d(value).astype(self.source_dtype, copy=COPY_IF_NEEDED)
 
         scaled = np.round(value * self.scale) + self.offset
 
@@ -423,7 +425,7 @@ def __call__(self, value):
         return result
 
     def inverse(self, value):
-        is_scalar = np.array(value, copy=False).shape == ()
+        is_scalar = np.array(value, copy=COPY_IF_NEEDED).shape == ()
         value = np.atleast_1d(value)
 
         result = (value.astype(self.source_dtype) - self.offset) / self.scale
@@ -530,7 +532,7 @@ def inverse(self, value):
 
         # astropy table columns somehow try to handle byte columns as strings
         # when iterating, this does not work here, convert to np.array
-        value = np.array(value, copy=False)
+        value = np.array(value, copy=COPY_IF_NEEDED)
         return np.array([v.decode("utf-8") for v in value])
 
     def get_meta(self, colname):

src/ctapipe/reco/impact.py

@@ -14,6 +14,7 @@
 
 from ctapipe.core import traits
 
+from ..compat import COPY_IF_NEEDED
 from ..containers import ReconstructedEnergyContainer, ReconstructedGeometryContainer
 from ..coordinates import (
     CameraFrame,
@@ -518,8 +519,8 @@ def get_likelihood(
         pix_x_rot, pix_y_rot = rotate_translate(
             self.pixel_y.data, self.pixel_x.data, source_y, source_x, -phi
         )
-        pix_x_rot = np.ma.array(pix_x_rot, mask=self.pixel_x.mask, copy=False)
-        pix_y_rot = np.ma.array(pix_y_rot, mask=self.pixel_y.mask, copy=False)
+        pix_x_rot = np.ma.array(pix_x_rot, mask=self.pixel_x.mask, copy=COPY_IF_NEEDED)
+        pix_y_rot = np.ma.array(pix_y_rot, mask=self.pixel_y.mask, copy=COPY_IF_NEEDED)
 
         # In the interpolator class we can gain speed advantages by using masked arrays
         # so we need to make sure here everything is masked

src/ctapipe/reco/sklearn.py

@@ -4,6 +4,7 @@
 import pathlib
 from abc import abstractmethod
 from collections import defaultdict
+from copy import deepcopy
 
 import astropy.units as u
 import joblib
@@ -57,6 +58,12 @@
 SUPPORTED_REGRESSORS = dict(all_estimators("regressor"))
 SUPPORTED_MODELS = {**SUPPORTED_CLASSIFIERS, **SUPPORTED_REGRESSORS}
 
+_invalid_geometry = ReconstructedGeometryContainer(
+    alt=u.Quantity(np.nan, unit=u.deg),
+    az=u.Quantity(np.nan, unit=u.deg),
+    is_valid=False,
+)
+
 
 class MLQualityQuery(QualityQuery):
     """Quality criteria for machine learning models with different defaults"""
@@ -758,20 +765,12 @@ def __call__(self, event: ArrayEventContainer) -> None:
                     disp_container = DispContainer(
                         parameter=u.Quantity(np.nan, self.unit),
                     )
-                    altaz_container = ReconstructedGeometryContainer(
-                        alt=u.Quantity(np.nan, u.deg, copy=False),
-                        az=u.Quantity(np.nan, u.deg, copy=False),
-                        is_valid=False,
-                    )
+                    altaz_container = deepcopy(_invalid_geometry)
             else:
                 disp_container = DispContainer(
                     parameter=u.Quantity(np.nan, self.unit),
                 )
-                altaz_container = ReconstructedGeometryContainer(
-                    alt=u.Quantity(np.nan, u.deg, copy=False),
-                    az=u.Quantity(np.nan, u.deg, copy=False),
-                    is_valid=False,
-                )
+                altaz_container = deepcopy(_invalid_geometry)
 
             disp_container.prefix = f"{self.prefix}_tel"
             altaz_container.prefix = f"{self.prefix}_tel"

src/ctapipe/reco/stereo_combination.py

@@ -10,6 +10,7 @@
 from ctapipe.core.traits import Bool, CaselessStrEnum, Unicode
 from ctapipe.reco.reconstructor import ReconstructionProperty
 
+from ..compat import COPY_IF_NEEDED
 from ..containers import (
     ArrayEventContainer,
     ParticleClassificationContainer,
@@ -186,8 +187,8 @@ def _combine_energy(self, event):
             valid = False
 
         event.dl2.stereo.energy[self.prefix] = ReconstructedEnergyContainer(
-            energy=u.Quantity(mean, u.TeV, copy=False),
-            energy_uncert=u.Quantity(std, u.TeV, copy=False),
+            energy=u.Quantity(mean, u.TeV, copy=COPY_IF_NEEDED),
+            energy_uncert=u.Quantity(std, u.TeV, copy=COPY_IF_NEEDED),
             telescopes=ids,
             is_valid=valid,
             prefix=self.prefix,
@@ -252,16 +253,16 @@ def _combine_altaz(self, event):
             valid = True
         else:
             mean_altaz = AltAz(
-                alt=u.Quantity(np.nan, u.deg, copy=False),
-                az=u.Quantity(np.nan, u.deg, copy=False),
+                alt=u.Quantity(np.nan, u.deg, copy=COPY_IF_NEEDED),
+                az=u.Quantity(np.nan, u.deg, copy=COPY_IF_NEEDED),
             )
             std = np.nan
             valid = False
 
         event.dl2.stereo.geometry[self.prefix] = ReconstructedGeometryContainer(
             alt=mean_altaz.alt,
             az=mean_altaz.az,
-            ang_distance_uncert=u.Quantity(np.rad2deg(std), u.deg, copy=False),
+            ang_distance_uncert=u.Quantity(np.rad2deg(std), u.deg, copy=COPY_IF_NEEDED),
             telescopes=ids,
             is_valid=valid,
             prefix=self.prefix,
@@ -358,11 +359,11 @@ def predict_table(self, mono_predictions: Table) -> Table:
                 std = np.full(n_array_events, np.nan)
 
             stereo_table[f"{self.prefix}_energy"] = u.Quantity(
-                stereo_energy, u.TeV, copy=False
+                stereo_energy, u.TeV, copy=COPY_IF_NEEDED
             )
 
             stereo_table[f"{self.prefix}_energy_uncert"] = u.Quantity(
-                std, u.TeV, copy=False
+                std, u.TeV, copy=COPY_IF_NEEDED
             )
             stereo_table[f"{self.prefix}_is_valid"] = np.isfinite(stereo_energy)
             stereo_table[f"{self.prefix}_goodness_of_fit"] = np.nan
@@ -397,16 +398,20 @@ def predict_table(self, mono_predictions: Table) -> Table:
                 std = np.sqrt(-2 * np.log(r))
             else:
                 mean_altaz = AltAz(
-                    alt=u.Quantity(np.full(n_array_events, np.nan), u.deg, copy=False),
-                    az=u.Quantity(np.full(n_array_events, np.nan), u.deg, copy=False),
+                    alt=u.Quantity(
+                        np.full(n_array_events, np.nan), u.deg, copy=COPY_IF_NEEDED
+                    ),
+                    az=u.Quantity(
+                        np.full(n_array_events, np.nan), u.deg, copy=COPY_IF_NEEDED
+                    ),
                 )
                 std = np.full(n_array_events, np.nan)
 
             stereo_table[f"{self.prefix}_alt"] = mean_altaz.alt.to(u.deg)
             stereo_table[f"{self.prefix}_az"] = mean_altaz.az.to(u.deg)
 
             stereo_table[f"{self.prefix}_ang_distance_uncert"] = u.Quantity(
-                np.rad2deg(std), u.deg, copy=False
+                np.rad2deg(std), u.deg, copy=COPY_IF_NEEDED
             )
 
             stereo_table[f"{self.prefix}_is_valid"] = np.logical_and(