diff --git a/.github/workflows/sonar.yml b/.github/workflows/sonar.yml index e31507b0662..0919c25e296 100644 --- a/.github/workflows/sonar.yml +++ b/.github/workflows/sonar.yml @@ -58,7 +58,7 @@ jobs: cat coverage/sonar_env - name: Sonarqube Scan - uses: SonarSource/sonarqube-scan-action@v2 + uses: SonarSource/sonarqube-scan-action@v3 env: GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} SONAR_TOKEN: ${{ secrets.SONAR_TOKEN }} diff --git a/.mailmap b/.mailmap index 37c9074b1e4..f7b247d65ee 100644 --- a/.mailmap +++ b/.mailmap @@ -93,6 +93,8 @@ Anno Knierim <32365583+aknierim@users.noreply.gith Vadym Voitsekhovskyi + + Tjark Miener Michael Punch diff --git a/AUTHORS b/AUTHORS index 4a17240000d..572bffc4fd1 100644 --- a/AUTHORS +++ b/AUTHORS @@ -11,10 +11,12 @@ Lukas Beiske Georg Schwefer Jonas Hackfeld Alison Mitchell +Tjark Miener Dominik Neise Michele Peresano Christoph Deil Justus Zorn +Christoph Toennis Anno Knierim Stefan Fröse Rune Michael Dominik @@ -24,12 +26,12 @@ Franca Cassol Clara Escanuela Nieves Thomas Vuillaume Kai Brügge +Mykhailo Dalchenko Felix Werner Satoshi Fukami Michele Mastropietro Jeremie DECOCK Abelardo Moralejo Olaizola -Mykhailo Dalchenko Tarek Hassan Wrijupan Bhattacharyya Francesco Visconti @@ -44,6 +46,7 @@ Tristan Carel Daniel Nieto Michael Punch Miguel Nievas +Arnau Aguasca-Cabot Christophe COSSOU Cyril Alispach Daniel Morcuende @@ -51,7 +54,6 @@ Julien Lefaucheur Konstantin Pfrang Moritz Hütten Thomas Armstrong -Tjark Miener Alice Donini <40267450+adonini@users.noreply.github.com> Andrés Baquero David Landriu diff --git a/CHANGES.rst b/CHANGES.rst index 74f97418edf..82d906de031 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -1,3 +1,60 @@ +ctapipe v0.22.0 (2024-09-12) +============================ + +API Changes +----------- + +- The ``PointingInterpolator`` was moved from ``ctapipe.io`` to ``ctapipe.monitoring``. [`#2615 `__] + + +Bug Fixes +--------- + +- Fix a redundant error message in ``Tool`` caused by normal ``SystemExit(0)`` [`#2575 `__] + +- Fix error message for non-existent config files. [`#2591 `__] + + +New Features +------------ + +- ctapipe is now compatible with numpy 2.0. [`#2580 `__] + Note: not all new behaviour of numpy 2.0 is followed, as the core dependency ``numba`` does not yet implement + all changes from numpy 2.0. See `the numba announcement for more detail `_. + +- Add lstchains image cleaning procedure including its pedestal cleaning method. [`#2541 `__] + +- A new ImageExtractor called ``VarianceExtractor`` was added + An Enum class was added to containers.py that is used in the metadata of the VarianceExtractor output [`#2543 `__] + +- Add API to extract the statistics from a sequence of images. [`#2554 `__] + +- The provenance system now records the reference metadata + of input and output files, if available. [`#2598 `__] + +- Add Interpolator class to generalize the PointingInterpolator in the monitoring collection. [`#2600 `__] + +- Add outlier detection components to identify faulty pixels. [`#2604 `__] + +- The ``ctapipe-merge`` tool now checks for duplicated input files and + raises an error in that case. + + The ``HDF5Merger`` class, and thus also the ``ctapipe-merge`` tool, + now checks for duplicated obs_ids during merging, to prevent + invalid output files. [`#2611 `__] + +- The ``Instrument.site`` metadata item now accepts any string, + not just a pre-defined list of sites. [`#2616 `__] + +Refactoring and Optimization +---------------------------- + +- Update exception handling in tools + + - Add a possibility to handle custom exception in ``Tool.run()`` + with the preservation of the exit code. [`#2594 `__] + + ctapipe v0.21.2 (2024-06-26) ============================ diff --git a/docs/_static/switcher.json b/docs/_static/switcher.json index f4eacc010e8..10a74c6763c 100644 --- a/docs/_static/switcher.json +++ b/docs/_static/switcher.json @@ -9,6 +9,11 @@ "version": "stable", "url": "https://ctapipe.readthedocs.io/en/stable/" }, + { + "name": "v0.22.0", + "version": "v0.22.0", + "url": "https://ctapipe.readthedocs.io/en/v0.22.0/" + }, { "name": "v0.21.2", "version": "v0.21.2", diff --git a/docs/api-reference/monitoring/index.rst b/docs/api-reference/monitoring/index.rst index 11ae54bd7e0..e34072bd6ce 100644 --- a/docs/api-reference/monitoring/index.rst +++ b/docs/api-reference/monitoring/index.rst @@ -18,10 +18,10 @@ Submodules .. toctree:: :maxdepth: 1 - :glob: aggregator calculator + interpolation outlier diff --git a/docs/api-reference/monitoring/interpolation.rst b/docs/api-reference/monitoring/interpolation.rst new file mode 100644 index 00000000000..7fab95f007f --- /dev/null +++ b/docs/api-reference/monitoring/interpolation.rst @@ -0,0 +1,11 @@ +.. _monitoring_interpolation: + +************* +Interpolation +************* + + +Reference/API +============= + +.. automodapi:: ctapipe.monitoring.interpolation diff --git a/docs/changes/2541.feature.rst b/docs/changes/2541.feature.rst deleted file mode 100644 index 30a27dbea21..00000000000 --- a/docs/changes/2541.feature.rst +++ /dev/null @@ -1 +0,0 @@ -Add lstchains image cleaning procedure including its pedestal cleaning method. diff --git a/docs/changes/2543.feature.rst b/docs/changes/2543.feature.rst deleted file mode 100644 index 8154498542d..00000000000 --- a/docs/changes/2543.feature.rst +++ /dev/null @@ -1,2 +0,0 @@ -A new ImageExtractor called ``VarianceExtractor`` was added -An Enum class was added to containers.py that is used in the metadata of the VarianceExtractor output diff --git a/docs/changes/2554.feature.rst b/docs/changes/2554.feature.rst deleted file mode 100644 index 2e6a6356b3a..00000000000 --- a/docs/changes/2554.feature.rst +++ /dev/null @@ -1 +0,0 @@ -Add API to extract the statistics from a sequence of images. diff --git a/docs/changes/2575.bugfix.rst b/docs/changes/2575.bugfix.rst deleted file mode 100644 index 6f6d0402165..00000000000 --- a/docs/changes/2575.bugfix.rst +++ /dev/null @@ -1 +0,0 @@ -Fix a redundant error message in ``Tool`` caused by normal ``SystemExit(0)`` diff --git a/docs/changes/2594.optimization.rst b/docs/changes/2594.optimization.rst deleted file mode 100644 index c1045f4d33c..00000000000 --- a/docs/changes/2594.optimization.rst +++ /dev/null @@ -1,4 +0,0 @@ -Update exception handling in tools - -- Add a possibility to handle custom exception in ``Tool.run()`` - with the preservation of the exit code. diff --git a/docs/changes/2598.feature.rst b/docs/changes/2598.feature.rst deleted file mode 100644 index 16f05623c72..00000000000 --- a/docs/changes/2598.feature.rst +++ /dev/null @@ -1,2 +0,0 @@ -The provenance system now records the reference metadata -of input and output files, if available. diff --git a/docs/changes/2602.bugfix.rst b/docs/changes/2602.bugfix.rst deleted file mode 100644 index 342758154ba..00000000000 --- a/docs/changes/2602.bugfix.rst +++ /dev/null @@ -1,3 +0,0 @@ -Fix ``NSBImageCleaner`` that tries to reach the tel attribute -from a MonitoringCameraContainer which results in an AttributeError -when using the ``ImageProcessor``. diff --git a/docs/changes/2604.feature.rst b/docs/changes/2604.feature.rst deleted file mode 100644 index 751bd59190a..00000000000 --- a/docs/changes/2604.feature.rst +++ /dev/null @@ -1 +0,0 @@ -Add outlier detection components to identify faulty pixels. diff --git a/docs/changes/template.rst b/docs/changes/template.rst index 7b4dfcb21b9..c9f1e06e24e 100644 --- a/docs/changes/template.rst +++ b/docs/changes/template.rst @@ -1,6 +1,6 @@ {% if render_title %} {% if versiondata.name %} -{{ versiondata.name }} {{ versiondata.version }} ({{ versiondata.date }}) +{{ versiondata.name | lower }} {{ versiondata.version }} ({{ versiondata.date }}) {{ top_underline * ((versiondata.name + versiondata.version + versiondata.date)|length + 4)}} {% else %} {{ versiondata.version }} ({{ versiondata.date }}) diff --git a/pyproject.toml b/pyproject.toml index 7a663d17c39..6c81e30568a 100644 --- a/pyproject.toml +++ b/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] diff --git a/sonar-project.properties b/sonar-project.properties index 12d5055665e..08e0dd7b9b0 100644 --- a/sonar-project.properties +++ b/sonar-project.properties @@ -5,4 +5,4 @@ sonar.python.version=3.10 # ignore examples for coverage and issues, these are sphinx-gallery notebook scripts # which aren't really supported by sonarqube and lead to a lot of false positives -sonar.exclusions=examples/** +sonar.exclusions=examples/**,coverage/coverage.xml 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/conftest.py b/src/ctapipe/conftest.py index 4d743e5d217..c6599dfe86d 100644 --- a/src/ctapipe/conftest.py +++ b/src/ctapipe/conftest.py @@ -592,6 +592,14 @@ def proton_train_clf(model_tmp_path, energy_regressor_path): ], raises=True, ) + + # modify obs_ids by adding a constant, this enables merging gamma and proton files + # which is used in the merge tool tests. + with tables.open_file(outpath, mode="r+") as f: + for table in f.walk_nodes("/", "Table"): + if "obs_id" in table.colnames: + obs_id = table.col("obs_id") + table.modify_column(colname="obs_id", column=obs_id + 1_000_000_000) return outpath 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/core/tool.py b/src/ctapipe/core/tool.py index 83fe9fcf28a..e484c6ef71d 100644 --- a/src/ctapipe/core/tool.py +++ b/src/ctapipe/core/tool.py @@ -150,15 +150,15 @@ def main(): trait=Path( exists=True, directory_ok=False, - help=( - "List of configuration files with parameters to load " - "in addition to command-line parameters. " - "The order listed is the order of precedence (later config parameters " - "overwrite earlier ones), however parameters specified on the " - "command line always have the highest precedence. " - "Config files may be in JSON, YAML, TOML, or Python format" - ), - ) + ), + help=( + "List of configuration files with parameters to load " + "in addition to command-line parameters. " + "The order listed is the order of precedence (later config parameters " + "overwrite earlier ones), however parameters specified on the " + "command line always have the highest precedence. " + "Config files may be in JSON, YAML, TOML, or Python format" + ), ).tag(config=True) log_config = Dict(default_value={}).tag(config=True) diff --git a/src/ctapipe/core/traits.py b/src/ctapipe/core/traits.py index cc14a8936bf..b787c3fad22 100644 --- a/src/ctapipe/core/traits.py +++ b/src/ctapipe/core/traits.py @@ -52,6 +52,17 @@ logger = logging.getLogger(__name__) +class PathError: + """Signal Non-Existence of a Path""" + + def __init__(self, path, reason): + self.path = path + self.reason = reason + + def __repr__(self): + return f"'{self.path}': {self.reason}" + + # Aliases Bool = traitlets.Bool Int = traitlets.Int @@ -111,9 +122,9 @@ def validate(self, obj, value): try: quantity = u.Quantity(value) except TypeError: - return self.error(obj, value) + self.error(obj, value) except ValueError: - return self.error(obj, value) + self.error(obj, value) if self.physical_type is not None: given_type = u.get_physical_type(quantity) @@ -136,7 +147,7 @@ def validate(self, obj, value): the_time.format = "iso" return the_time except ValueError: - return self.error(obj, value) + self.error(obj, value) def info(self): info = "an ISO8601 datestring or Time instance" @@ -204,19 +215,19 @@ def validate(self, obj, value): self.error(obj, value) if not isinstance(value, str | pathlib.Path): - return self.error(obj, value) + self.error(obj, value) # expand any environment variables in the path: value = os.path.expandvars(value) if isinstance(value, str): if value == "": - return self.error(obj, value) + self.error(obj, value) try: url = urlparse(value) except ValueError: - return self.error(obj, value) + self.error(obj, value) if url.scheme in ("http", "https"): # here to avoid circular import, since every module imports @@ -233,23 +244,18 @@ def validate(self, obj, value): elif url.scheme in ("", "file"): value = pathlib.Path(url.netloc, url.path) else: - return self.error(obj, value) + self.error(obj, value) value = value.absolute() exists = value.exists() if self.exists is not None: if exists != self.exists: - raise TraitError( - 'Path "{}" {} exist'.format( - value, "does not" if self.exists else "must not" - ) - ) + raise TraitError(PathError(value, "does not exist"), self.info(), self) if exists: if not self.directory_ok and value.is_dir(): - raise TraitError(f'Path "{value}" must not be a directory') + raise TraitError(PathError(value, "is a directory"), self.info(), self) if not self.file_ok and value.is_file(): - raise TraitError(f'Path "{value}" must not be a file') - + raise TraitError(PathError(value, "is a file"), self.info(), self) return value diff --git a/src/ctapipe/image/extractor.py b/src/ctapipe/image/extractor.py index 88e7e683e9b..a8a6a71323c 100644 --- a/src/ctapipe/image/extractor.py +++ b/src/ctapipe/image/extractor.py @@ -740,7 +740,7 @@ def __call__( class NeighborPeakWindowSum(ImageExtractor): """ Extractor which sums in a window about the - peak defined by the wavefroms in neighboring pixels. + peak defined by the waveforms in neighboring pixels. """ window_width = IntTelescopeParameter( @@ -808,7 +808,7 @@ def __call__( class BaselineSubtractedNeighborPeakWindowSum(NeighborPeakWindowSum): """ Extractor that first subtracts the baseline before summing in a - window about the peak defined by the wavefroms in neighboring pixels. + window about the peak defined by the waveforms in neighboring pixels. """ baseline_start = Int(0, help="Start sample for baseline estimation").tag( @@ -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/camera/geometry.py b/src/ctapipe/instrument/camera/geometry.py index ab8aa84e94d..d1f861d7549 100644 --- a/src/ctapipe/instrument/camera/geometry.py +++ b/src/ctapipe/instrument/camera/geometry.py @@ -76,7 +76,7 @@ def from_string(cls, name): class CameraGeometry: """`CameraGeometry` is a class that stores information about a - Cherenkov Camera that us useful for imaging algorithms and + Cherenkov Camera that is useful for imaging algorithms and displays. It contains lists of pixel positions, areas, pixel shapes, as well as a neighbor (adjacency) list and matrix for each pixel. In general the neighbor_matrix attribute should be used in any algorithm 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/__init__.py b/src/ctapipe/io/__init__.py index afe96f39430..cfb5fc5501e 100644 --- a/src/ctapipe/io/__init__.py +++ b/src/ctapipe/io/__init__.py @@ -18,7 +18,6 @@ from .datawriter import DATA_MODEL_VERSION, DataWriter - __all__ = [ "HDF5TableWriter", "HDF5TableReader", diff --git a/src/ctapipe/io/hdf5eventsource.py b/src/ctapipe/io/hdf5eventsource.py index afbe2153d2c..0cb47d1aee8 100644 --- a/src/ctapipe/io/hdf5eventsource.py +++ b/src/ctapipe/io/hdf5eventsource.py @@ -7,9 +7,7 @@ from astropy.table import QTable from astropy.utils.decorators import lazyproperty -from ctapipe.atmosphere import AtmosphereDensityProfile -from ctapipe.instrument.optics import FocalLengthKind - +from ..atmosphere import AtmosphereDensityProfile from ..containers import ( ArrayEventContainer, CameraHillasParametersContainer, @@ -48,11 +46,12 @@ from ..core import Container, Field from ..core.traits import UseEnum from ..instrument import SubarrayDescription +from ..instrument.optics import FocalLengthKind +from ..monitoring.interpolation import PointingInterpolator from .astropy_helpers import read_table from .datalevels import DataLevel from .eventsource import EventSource from .hdf5tableio import HDF5TableReader -from .pointing import PointingInterpolator from .tableloader import DL2_SUBARRAY_GROUP, DL2_TELESCOPE_GROUP, POINTING_GROUP __all__ = ["HDF5EventSource"] diff --git a/src/ctapipe/io/hdf5merger.py b/src/ctapipe/io/hdf5merger.py index 95d0a8b6909..b6d59ae7364 100644 --- a/src/ctapipe/io/hdf5merger.py +++ b/src/ctapipe/io/hdf5merger.py @@ -188,6 +188,8 @@ def __init__(self, output_path=None, **kwargs): self.data_model_version = None self.subarray = None self.meta = None + self._merged_obs_ids = set() + # output file existed, so read subarray and data model version to make sure # any file given matches what we already have if appending: @@ -202,6 +204,9 @@ def __init__(self, output_path=None, **kwargs): ) self.required_nodes = _get_required_nodes(self.h5file) + # this will update _merged_obs_ids from existing input file + self._check_obs_ids(self.h5file) + def __call__(self, other: str | Path | tables.File): """ Append file ``other`` to the output file @@ -267,7 +272,32 @@ def _check_can_merge(self, other): f"Required node {node_path} not found in {other.filename}" ) + def _check_obs_ids(self, other): + keys = [ + "/configuration/observation/observation_block", + "/dl1/event/subarray/trigger", + ] + + for key in keys: + if key in other.root: + obs_ids = other.root[key].col("obs_id") + break + else: + raise CannotMerge( + f"Input file {other.filename} is missing keys required to" + f" check for duplicated obs_ids. Tried: {keys}" + ) + + duplicated = self._merged_obs_ids.intersection(obs_ids) + if len(duplicated) > 0: + msg = f"Input file {other.filename} contains obs_ids already included in output file: {duplicated}" + raise CannotMerge(msg) + + self._merged_obs_ids.update(obs_ids) + def _append(self, other): + self._check_obs_ids(other) + # Configuration self._append_subarray(other) diff --git a/src/ctapipe/io/metadata.py b/src/ctapipe/io/metadata.py index 2aef43d111e..b7211bc0265 100644 --- a/src/ctapipe/io/metadata.py +++ b/src/ctapipe/io/metadata.py @@ -176,23 +176,14 @@ def default_time(self): class Instrument(Configurable): """Instrumental Context""" - site = Enum( - [ - "CTA-North", - "CTA-South", - "SDMC", - "HQ", - "MAGIC", - "HESS", - "VERITAS", - "FACT", - "Whipple", - "Other", - ], - "Other", - help="Which site of CTA (or external telescope) " - "this instrument is associated with", + site = Unicode( + default_value="Other", + help=( + "Which site of CTAO (or external telescope)" + " this instrument is associated with" + ), ).tag(config=True) + class_ = Enum( [ "Array", diff --git a/src/ctapipe/io/pointing.py b/src/ctapipe/io/pointing.py deleted file mode 100644 index 412c0510a98..00000000000 --- a/src/ctapipe/io/pointing.py +++ /dev/null @@ -1,137 +0,0 @@ -from typing import Any - -import astropy.units as u -import numpy as np -import tables -from astropy.time import Time -from scipy.interpolate import interp1d - -from ctapipe.core import Component, traits - -from .astropy_helpers import read_table - - -class PointingInterpolator(Component): - """ - Interpolate pointing from a monitoring table to a given timestamp. - - Parameters - ---------- - h5file : None | tables.File - A open hdf5 file with read access. - The monitoring table is expected to be stored in that file at - ``/dl0/monitoring/telescope/pointing/tel_{tel_id:03d}`` - - If not given, monitoring tables can be added via `PointingInterpolator.add_table`. - """ - - bounds_error = traits.Bool( - default_value=True, - help="If true, raises an exception when trying to extrapolate out of the given table", - ).tag(config=True) - - extrapolate = traits.Bool( - help="If bounds_error is False, this flag will specify whether values outside" - "the available values are filled with nan (False) or extrapolated (True).", - default_value=False, - ).tag(config=True) - - def __init__(self, h5file=None, **kwargs): - super().__init__(**kwargs) - - if h5file is not None and not isinstance(h5file, tables.File): - raise TypeError("h5file must be a tables.File") - self.h5file = h5file - - self.interp_options: dict[str, Any] = dict(assume_sorted=True, copy=False) - if self.bounds_error: - self.interp_options["bounds_error"] = True - elif self.extrapolate: - self.interp_options["bounds_error"] = False - self.interp_options["fill_value"] = "extrapolate" - else: - self.interp_options["bounds_error"] = False - self.interp_options["fill_value"] = np.nan - - self._alt_interpolators = {} - self._az_interpolators = {} - - def add_table(self, tel_id, pointing_table): - """ - Add a table to this interpolator - - Parameters - ---------- - tel_id : int - Telescope id - pointing_table : astropy.table.Table - Table of pointing values, expected columns - are ``time`` as ``Time`` column, ``azimuth`` and ``altitude`` - as quantity columns. - """ - missing = {"time", "azimuth", "altitude"} - set(pointing_table.colnames) - if len(missing) > 0: - raise ValueError(f"Table is missing required column(s): {missing}") - - if not isinstance(pointing_table["time"], Time): - raise TypeError("'time' column of pointing table must be astropy.time.Time") - - for col in ("azimuth", "altitude"): - unit = pointing_table[col].unit - if unit is None or not u.rad.is_equivalent(unit): - raise ValueError(f"{col} must have units compatible with 'rad'") - - # sort first, so it's not done twice for each interpolator - pointing_table.sort("time") - # interpolate in mjd TAI. Float64 mjd is precise enough for pointing - # and TAI is contiguous, so no issues with leap seconds. - mjd = pointing_table["time"].tai.mjd - - az = pointing_table["azimuth"].quantity.to_value(u.rad) - # prepare azimuth for interpolation by "unwrapping": i.e. turning - # [359, 1] into [359, 361]. This assumes that if we get values like - # [359, 1] the telescope moved 2 degrees through 0, not 358 degrees - # the other way around. This should be true for all telescopes given - # the sampling speed of pointing values and their maximum movement speed. - # No telescope can turn more than 180° in 2 seconds. - az = np.unwrap(az) - alt = pointing_table["altitude"].quantity.to_value(u.rad) - - self._az_interpolators[tel_id] = interp1d(mjd, az, **self.interp_options) - self._alt_interpolators[tel_id] = interp1d(mjd, alt, **self.interp_options) - - def _read_pointing_table(self, tel_id): - pointing_table = read_table( - self.h5file, - f"/dl0/monitoring/telescope/pointing/tel_{tel_id:03d}", - ) - self.add_table(tel_id, pointing_table) - - def __call__(self, tel_id, time): - """ - Interpolate alt/az for given time and tel_id. - - Parameters - ---------- - tel_id : int - telescope id - time : astropy.time.Time - time for which to interpolate the pointing - - Returns - ------- - altitude : astropy.units.Quantity[deg] - interpolated altitude angle - azimuth : astropy.units.Quantity[deg] - interpolated azimuth angle - """ - if tel_id not in self._az_interpolators: - if self.h5file is not None: - self._read_pointing_table(tel_id) - else: - raise KeyError(f"No pointing table available for tel_id {tel_id}") - - mjd = time.tai.mjd - az = u.Quantity(self._az_interpolators[tel_id](mjd), u.rad, copy=False) - alt = u.Quantity(self._alt_interpolators[tel_id](mjd), u.rad, copy=False) - return alt, az 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/io/tableloader.py b/src/ctapipe/io/tableloader.py index 36a954fb7de..2b10ad96506 100644 --- a/src/ctapipe/io/tableloader.py +++ b/src/ctapipe/io/tableloader.py @@ -13,8 +13,8 @@ from ..core import Component, Provenance, traits from ..instrument import FocalLengthKind, SubarrayDescription +from ..monitoring.interpolation import PointingInterpolator from .astropy_helpers import join_allow_empty, read_table -from .pointing import PointingInterpolator __all__ = ["TableLoader"] diff --git a/src/ctapipe/io/tests/test_merge.py b/src/ctapipe/io/tests/test_merge.py index ec37b6b5dda..6376da6566a 100644 --- a/src/ctapipe/io/tests/test_merge.py +++ b/src/ctapipe/io/tests/test_merge.py @@ -68,7 +68,7 @@ def test_simple(tmp_path, gamma_train_clf, proton_train_clf): merger(proton_train_clf) subarray = SubarrayDescription.from_hdf(gamma_train_clf) - assert subarray == SubarrayDescription.from_hdf(output), "Subarays do not match" + assert subarray == SubarrayDescription.from_hdf(output), "Subarrays do not match" tel_groups = [ "/dl1/event/telescope/parameters", @@ -164,3 +164,28 @@ def test_muon(tmp_path, dl1_muon_output_file): n_input = len(input_table) assert len(table) == n_input assert_table_equal(table, input_table) + + +def test_duplicated_obs_ids(tmp_path, dl2_shower_geometry_file): + from ctapipe.io.hdf5merger import CannotMerge, HDF5Merger + + output = tmp_path / "invalid.dl1.h5" + + # check for fresh file + with HDF5Merger(output) as merger: + merger(dl2_shower_geometry_file) + + with pytest.raises( + CannotMerge, match="Input file .* contains obs_ids already included" + ): + merger(dl2_shower_geometry_file) + + # check for appending + with HDF5Merger(output, overwrite=True) as merger: + merger(dl2_shower_geometry_file) + + with HDF5Merger(output, append=True) as merger: + with pytest.raises( + CannotMerge, match="Input file .* contains obs_ids already included" + ): + merger(dl2_shower_geometry_file) diff --git a/src/ctapipe/monitoring/__init__.py b/src/ctapipe/monitoring/__init__.py new file mode 100644 index 00000000000..cb102f768cf --- /dev/null +++ b/src/ctapipe/monitoring/__init__.py @@ -0,0 +1,23 @@ +""" +Module for handling monitoring data. +""" +from .aggregator import PlainAggregator, SigmaClippingAggregator, StatisticsAggregator +from .interpolation import Interpolator, PointingInterpolator +from .outlier import ( + MedianOutlierDetector, + OutlierDetector, + RangeOutlierDetector, + StdOutlierDetector, +) + +__all__ = [ + "PlainAggregator", + "SigmaClippingAggregator", + "StatisticsAggregator", + "OutlierDetector", + "RangeOutlierDetector", + "MedianOutlierDetector", + "StdOutlierDetector", + "Interpolator", + "PointingInterpolator", +] diff --git a/src/ctapipe/monitoring/aggregator.py b/src/ctapipe/monitoring/aggregator.py index a6c8be44c5f..8102b21e855 100644 --- a/src/ctapipe/monitoring/aggregator.py +++ b/src/ctapipe/monitoring/aggregator.py @@ -25,6 +25,12 @@ from ctapipe.core import TelescopeComponent from ctapipe.core.traits import Int +__all__ = [ + "StatisticsAggregator", + "PlainAggregator", + "SigmaClippingAggregator", +] + class StatisticsAggregator(TelescopeComponent): """ diff --git a/src/ctapipe/monitoring/interpolation.py b/src/ctapipe/monitoring/interpolation.py new file mode 100644 index 00000000000..84064cbc1a3 --- /dev/null +++ b/src/ctapipe/monitoring/interpolation.py @@ -0,0 +1,188 @@ +from abc import ABCMeta, abstractmethod +from typing import Any + +import astropy.units as u +import numpy as np +import tables +from astropy.time import Time +from scipy.interpolate import interp1d + +from ctapipe.core import Component, traits + +__all__ = [ + "Interpolator", + "PointingInterpolator", +] + + +class Interpolator(Component, metaclass=ABCMeta): + """ + Interpolator parent class. + + Parameters + ---------- + h5file : None | tables.File + A open hdf5 file with read access. + """ + + bounds_error = traits.Bool( + default_value=True, + help="If true, raises an exception when trying to extrapolate out of the given table", + ).tag(config=True) + + extrapolate = traits.Bool( + help="If bounds_error is False, this flag will specify whether values outside" + "the available values are filled with nan (False) or extrapolated (True).", + default_value=False, + ).tag(config=True) + + telescope_data_group = None + required_columns = set() + expected_units = {} + + def __init__(self, h5file=None, **kwargs): + super().__init__(**kwargs) + + if h5file is not None and not isinstance(h5file, tables.File): + raise TypeError("h5file must be a tables.File") + self.h5file = h5file + + self.interp_options: dict[str, Any] = dict(assume_sorted=True, copy=False) + if self.bounds_error: + self.interp_options["bounds_error"] = True + elif self.extrapolate: + self.interp_options["bounds_error"] = False + self.interp_options["fill_value"] = "extrapolate" + else: + self.interp_options["bounds_error"] = False + self.interp_options["fill_value"] = np.nan + + self._interpolators = {} + + @abstractmethod + def add_table(self, tel_id, input_table): + """ + Add a table to this interpolator + This method reads input tables and creates instances of the needed interpolators + to be added to _interpolators. The first index of _interpolators needs to be + tel_id, the second needs to be the name of the parameter that is to be interpolated + + Parameters + ---------- + tel_id : int + Telescope id + input_table : astropy.table.Table + Table of pointing values, expected columns + are always ``time`` as ``Time`` column and + other columns for the data that is to be interpolated + """ + + pass + + def _check_tables(self, input_table): + missing = self.required_columns - set(input_table.colnames) + if len(missing) > 0: + raise ValueError(f"Table is missing required column(s): {missing}") + for col in self.expected_units: + unit = input_table[col].unit + if unit is None: + if self.expected_units[col] is not None: + raise ValueError( + f"{col} must have units compatible with '{self.expected_units[col].name}'" + ) + elif not self.expected_units[col].is_equivalent(unit): + if self.expected_units[col] is None: + raise ValueError(f"{col} must have units compatible with 'None'") + else: + raise ValueError( + f"{col} must have units compatible with '{self.expected_units[col].name}'" + ) + + def _check_interpolators(self, tel_id): + if tel_id not in self._interpolators: + if self.h5file is not None: + self._read_parameter_table(tel_id) # might need to be removed + else: + raise KeyError(f"No table available for tel_id {tel_id}") + + def _read_parameter_table(self, tel_id): + # prevent circular import between io and monitoring + from ..io import read_table + + input_table = read_table( + self.h5file, + f"{self.telescope_data_group}/tel_{tel_id:03d}", + ) + self.add_table(tel_id, input_table) + + +class PointingInterpolator(Interpolator): + """ + Interpolator for pointing and pointing correction data + """ + + telescope_data_group = "/dl0/monitoring/telescope/pointing" + required_columns = frozenset(["time", "azimuth", "altitude"]) + expected_units = {"azimuth": u.rad, "altitude": u.rad} + + def __call__(self, tel_id, time): + """ + Interpolate alt/az for given time and tel_id. + + Parameters + ---------- + tel_id : int + telescope id + time : astropy.time.Time + time for which to interpolate the pointing + + Returns + ------- + altitude : astropy.units.Quantity[deg] + interpolated altitude angle + azimuth : astropy.units.Quantity[deg] + interpolated azimuth angle + """ + + self._check_interpolators(tel_id) + + mjd = time.tai.mjd + az = u.Quantity(self._interpolators[tel_id]["az"](mjd), u.rad, copy=False) + alt = u.Quantity(self._interpolators[tel_id]["alt"](mjd), u.rad, copy=False) + return alt, az + + def add_table(self, tel_id, input_table): + """ + Add a table to this interpolator + + Parameters + ---------- + tel_id : int + Telescope id + input_table : astropy.table.Table + Table of pointing values, expected columns + are ``time`` as ``Time`` column, ``azimuth`` and ``altitude`` + as quantity columns for pointing and pointing correction data. + """ + + self._check_tables(input_table) + + if not isinstance(input_table["time"], Time): + raise TypeError("'time' column of pointing table must be astropy.time.Time") + + input_table = input_table.copy() + input_table.sort("time") + + az = input_table["azimuth"].quantity.to_value(u.rad) + # prepare azimuth for interpolation by "unwrapping": i.e. turning + # [359, 1] into [359, 361]. This assumes that if we get values like + # [359, 1] the telescope moved 2 degrees through 0, not 358 degrees + # the other way around. This should be true for all telescopes given + # the sampling speed of pointing values and their maximum movement speed. + # No telescope can turn more than 180° in 2 seconds. + az = np.unwrap(az) + alt = input_table["altitude"].quantity.to_value(u.rad) + mjd = input_table["time"].tai.mjd + self._interpolators[tel_id] = {} + self._interpolators[tel_id]["az"] = interp1d(mjd, az, **self.interp_options) + self._interpolators[tel_id]["alt"] = interp1d(mjd, alt, **self.interp_options) diff --git a/src/ctapipe/monitoring/outlier.py b/src/ctapipe/monitoring/outlier.py index 93b67bfd418..b2b84711cde 100644 --- a/src/ctapipe/monitoring/outlier.py +++ b/src/ctapipe/monitoring/outlier.py @@ -1,6 +1,12 @@ """ Outlier detection algorithms to identify faulty pixels """ +from abc import abstractmethod + +import numpy as np + +from ctapipe.core import TelescopeComponent +from ctapipe.core.traits import Float, List __all__ = [ "OutlierDetector", @@ -9,13 +15,6 @@ "StdOutlierDetector", ] -from abc import abstractmethod - -import numpy as np - -from ctapipe.core import TelescopeComponent -from ctapipe.core.traits import Float, List - class OutlierDetector(TelescopeComponent): """ diff --git a/src/ctapipe/io/tests/test_pointing.py b/src/ctapipe/monitoring/tests/test_interpolator.py similarity index 65% rename from src/ctapipe/io/tests/test_pointing.py rename to src/ctapipe/monitoring/tests/test_interpolator.py index 10913053c18..782aeae7435 100644 --- a/src/ctapipe/io/tests/test_pointing.py +++ b/src/ctapipe/monitoring/tests/test_interpolator.py @@ -5,37 +5,13 @@ from astropy.table import Table from astropy.time import Time +from ctapipe.monitoring.interpolation import PointingInterpolator -def test_simple(): - """Test pointing interpolation""" - from ctapipe.io.pointing import PointingInterpolator - - t0 = Time("2022-01-01T00:00:00") - - table = Table( - { - "time": t0 + np.arange(0.0, 10.1, 2.0) * u.s, - "azimuth": np.linspace(0.0, 10.0, 6) * u.deg, - "altitude": np.linspace(70.0, 60.0, 6) * u.deg, - }, - ) - - interpolator = PointingInterpolator() - interpolator.add_table(1, table) - - alt, az = interpolator(tel_id=1, time=t0 + 1 * u.s) - assert u.isclose(alt, 69 * u.deg) - assert u.isclose(az, 1 * u.deg) - - with pytest.raises(KeyError): - interpolator(tel_id=2, time=t0 + 1 * u.s) +t0 = Time("2022-01-01T00:00:00") def test_azimuth_switchover(): """Test pointing interpolation""" - from ctapipe.io.pointing import PointingInterpolator - - t0 = Time("2022-01-01T00:00:00") table = Table( { @@ -55,7 +31,6 @@ def test_azimuth_switchover(): def test_invalid_input(): """Test invalid pointing tables raise nice errors""" - from ctapipe.io.pointing import PointingInterpolator wrong_time = Table( { @@ -91,10 +66,8 @@ def test_invalid_input(): def test_hdf5(tmp_path): + """Test writing interpolated data to file""" from ctapipe.io import write_table - from ctapipe.io.pointing import PointingInterpolator - - t0 = Time("2022-01-01T00:00:00") table = Table( { @@ -115,11 +88,8 @@ def test_hdf5(tmp_path): def test_bounds(): """Test invalid pointing tables raise nice errors""" - from ctapipe.io.pointing import PointingInterpolator - t0 = Time("2022-01-01T00:00:00") - - table = Table( + table_pointing = Table( { "time": t0 + np.arange(0.0, 10.1, 2.0) * u.s, "azimuth": np.linspace(0.0, 10.0, 6) * u.deg, @@ -127,28 +97,36 @@ def test_bounds(): }, ) - interpolator = PointingInterpolator() - interpolator.add_table(1, table) + interpolator_pointing = PointingInterpolator() + interpolator_pointing.add_table(1, table_pointing) + error_message = "below the interpolation range" - with pytest.raises(ValueError, match="below the interpolation range"): - interpolator(tel_id=1, time=t0 - 0.1 * u.s) + with pytest.raises(ValueError, match=error_message): + interpolator_pointing(tel_id=1, time=t0 - 0.1 * u.s) with pytest.raises(ValueError, match="above the interpolation range"): - interpolator(tel_id=1, time=t0 + 10.2 * u.s) + interpolator_pointing(tel_id=1, time=t0 + 10.2 * u.s) - interpolator = PointingInterpolator(bounds_error=False) - interpolator.add_table(1, table) + alt, az = interpolator_pointing(tel_id=1, time=t0 + 1 * u.s) + assert u.isclose(alt, 69 * u.deg) + assert u.isclose(az, 1 * u.deg) + + with pytest.raises(KeyError): + interpolator_pointing(tel_id=2, time=t0 + 1 * u.s) + + interpolator_pointing = PointingInterpolator(bounds_error=False) + interpolator_pointing.add_table(1, table_pointing) for dt in (-0.1, 10.1) * u.s: - alt, az = interpolator(tel_id=1, time=t0 + dt) + alt, az = interpolator_pointing(tel_id=1, time=t0 + dt) assert np.isnan(alt.value) assert np.isnan(az.value) - interpolator = PointingInterpolator(bounds_error=False, extrapolate=True) - interpolator.add_table(1, table) - alt, az = interpolator(tel_id=1, time=t0 - 1 * u.s) + interpolator_pointing = PointingInterpolator(bounds_error=False, extrapolate=True) + interpolator_pointing.add_table(1, table_pointing) + alt, az = interpolator_pointing(tel_id=1, time=t0 - 1 * u.s) assert u.isclose(alt, 71 * u.deg) assert u.isclose(az, -1 * u.deg) - alt, az = interpolator(tel_id=1, time=t0 + 11 * u.s) + alt, az = interpolator_pointing(tel_id=1, time=t0 + 11 * u.s) assert u.isclose(alt, 59 * u.deg) assert u.isclose(az, 11 * u.deg) 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( diff --git a/src/ctapipe/tools/merge.py b/src/ctapipe/tools/merge.py index 3665d504cd4..11285ad9caf 100644 --- a/src/ctapipe/tools/merge.py +++ b/src/ctapipe/tools/merge.py @@ -3,6 +3,7 @@ """ import sys from argparse import ArgumentParser +from collections import Counter from pathlib import Path from tqdm.auto import tqdm @@ -161,6 +162,13 @@ def setup(self): ) sys.exit(1) + counts = Counter(self.input_files) + duplicated = [p for p, c in counts.items() if c > 1] + if len(duplicated) > 0: + raise ToolConfigurationError( + f"Same file given multiple times. Duplicated files are: {duplicated}" + ) + self.merger = self.enter_context(HDF5Merger(parent=self)) if self.merger.output_path in self.input_files: raise ToolConfigurationError( diff --git a/src/ctapipe/tools/tests/test_merge.py b/src/ctapipe/tools/tests/test_merge.py index 09b3bca217e..0879d27f1df 100644 --- a/src/ctapipe/tools/tests/test_merge.py +++ b/src/ctapipe/tools/tests/test_merge.py @@ -5,11 +5,12 @@ from pathlib import Path import numpy as np +import pytest import tables from astropy.table import vstack from astropy.utils.diff import report_diff_values -from ctapipe.core import run_tool +from ctapipe.core import ToolConfigurationError, run_tool from ctapipe.io import TableLoader from ctapipe.io.astropy_helpers import read_table from ctapipe.io.tests.test_astropy_helpers import assert_table_equal @@ -176,7 +177,6 @@ def test_muon(tmp_path, dl1_muon_output_file): argv=[ f"--output={output}", str(dl1_muon_output_file), - str(dl1_muon_output_file), ], raises=True, ) @@ -185,6 +185,24 @@ def test_muon(tmp_path, dl1_muon_output_file): input_table = read_table(dl1_muon_output_file, "/dl1/event/telescope/muon/tel_001") n_input = len(input_table) - assert len(table) == 2 * n_input - assert_table_equal(table[:n_input], input_table) - assert_table_equal(table[n_input:], input_table) + assert len(table) == n_input + assert_table_equal(table, input_table) + + +def test_duplicated(tmp_path, dl1_file, dl1_proton_file): + from ctapipe.tools.merge import MergeTool + + output = tmp_path / "invalid.dl1.h5" + with pytest.raises(ToolConfigurationError, match="Same file given multiple times"): + run_tool( + MergeTool(), + argv=[ + str(dl1_file), + str(dl1_proton_file), + str(dl1_file), + f"--output={output}", + "--overwrite", + ], + cwd=tmp_path, + raises=True, + )