From 89df66dbd9e04835fd1623b3be1e8e110e5f9f45 Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Wed, 4 Sep 2024 17:58:42 +0200 Subject: [PATCH] Add compatibility with numpy 2.0 --- pyproject.toml | 1 + src/ctapipe/atmosphere.py | 10 +++++++--- src/ctapipe/compat.py | 11 +++++++++++ src/ctapipe/coordinates/camera_frame.py | 9 +++++++-- src/ctapipe/image/extractor.py | 2 +- src/ctapipe/instrument/subarray.py | 9 +++++---- src/ctapipe/io/simteleventsource.py | 17 +++++++++++------ src/ctapipe/io/tableio.py | 14 ++++++++------ src/ctapipe/reco/impact.py | 5 +++-- src/ctapipe/reco/sklearn.py | 19 +++++++++---------- src/ctapipe/reco/stereo_combination.py | 25 +++++++++++++++---------- 11 files changed, 78 insertions(+), 44 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 04b8de8115c..6c81e30568a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -48,6 +48,7 @@ dependencies = [ "tqdm >=4.32", "traitlets ~=5.6", "zstandard", + "packaging", ] [project.optional-dependencies] diff --git a/src/ctapipe/atmosphere.py b/src/ctapipe/atmosphere.py index 0fcd865054c..62f25bc09ba 100644 --- a/src/ctapipe/atmosphere.py +++ b/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 ( diff --git a/src/ctapipe/compat.py b/src/ctapipe/compat.py index e9efad38758..34cb4e41bbf 100644 --- a/src/ctapipe/compat.py +++ b/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 diff --git a/src/ctapipe/coordinates/camera_frame.py b/src/ctapipe/coordinates/camera_frame.py index 766a3b41678..49771a78132 100644 --- a/src/ctapipe/coordinates/camera_frame.py +++ b/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) diff --git a/src/ctapipe/image/extractor.py b/src/ctapipe/image/extractor.py index a267a860f4b..a8a6a71323c 100644 --- a/src/ctapipe/image/extractor.py +++ b/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 diff --git a/src/ctapipe/instrument/subarray.py b/src/ctapipe/instrument/subarray.py index f4797aa7c1a..e69e34f77f3 100644 --- a/src/ctapipe/instrument/subarray.py +++ b/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) }, diff --git a/src/ctapipe/io/simteleventsource.py b/src/ctapipe/io/simteleventsource.py index cf89b11b461..175914fbbe9 100644 --- a/src/ctapipe/io/simteleventsource.py +++ b/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) diff --git a/src/ctapipe/io/tableio.py b/src/ctapipe/io/tableio.py index 9db162d5495..442cdaabcd1 100644 --- a/src/ctapipe/io/tableio.py +++ b/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): diff --git a/src/ctapipe/reco/impact.py b/src/ctapipe/reco/impact.py index b35c0af9a95..e837e8a8611 100644 --- a/src/ctapipe/reco/impact.py +++ b/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 diff --git a/src/ctapipe/reco/sklearn.py b/src/ctapipe/reco/sklearn.py index 19655dc1a0c..835339f3d08 100644 --- a/src/ctapipe/reco/sklearn.py +++ b/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" diff --git a/src/ctapipe/reco/stereo_combination.py b/src/ctapipe/reco/stereo_combination.py index 7c48a971fe5..19c561d89b2 100644 --- a/src/ctapipe/reco/stereo_combination.py +++ b/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,8 +253,8 @@ 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 @@ -261,7 +262,7 @@ def _combine_altaz(self, event): 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,8 +398,12 @@ 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) @@ -406,7 +411,7 @@ def predict_table(self, mono_predictions: Table) -> Table: 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(