From 829abc531b1435eb97770dadf78e228698b6b5a6 Mon Sep 17 00:00:00 2001 From: Lukas Beiske Date: Fri, 22 Nov 2024 17:17:35 +0100 Subject: [PATCH 1/8] Fix uncertainty calculation in StereoMeanCombiner and add vectorization module --- docs/api-reference/vectorization/index.rst | 20 +++ docs/changes/2658.bugfix.rst | 2 + src/ctapipe/reco/stereo_combination.py | 117 ++++++-------- .../reco/tests/test_stereo_combination.py | 10 ++ src/ctapipe/tools/apply_models.py | 34 ++-- src/ctapipe/vectorization/__init__.py | 6 + src/ctapipe/vectorization/aggregate.py | 145 ++++++++++++++++++ .../vectorization/tests/test_aggregate.py | 50 ++++++ 8 files changed, 300 insertions(+), 84 deletions(-) create mode 100644 docs/api-reference/vectorization/index.rst create mode 100644 docs/changes/2658.bugfix.rst create mode 100644 src/ctapipe/vectorization/__init__.py create mode 100644 src/ctapipe/vectorization/aggregate.py create mode 100644 src/ctapipe/vectorization/tests/test_aggregate.py diff --git a/docs/api-reference/vectorization/index.rst b/docs/api-reference/vectorization/index.rst new file mode 100644 index 00000000000..640a28b2eb0 --- /dev/null +++ b/docs/api-reference/vectorization/index.rst @@ -0,0 +1,20 @@ +.. _vectorization: + +******************************************************************************* +Helper Functions for Vectorized Operations on Tables (`~ctapipe.vectorization`) +******************************************************************************* + +.. currentmodule:: ctapipe.vectorization + +This module contains helper functions for performing vectorized operations +on tables, as used in e.g. :ref:`reco_stereo_combination`. + + +Reference/API +============= + +.. automodapi:: ctapipe.vectorization + :no-inheritance-diagram: + +.. automodapi:: ctapipe.vectorization.aggregate + :no-inheritance-diagram: diff --git a/docs/changes/2658.bugfix.rst b/docs/changes/2658.bugfix.rst new file mode 100644 index 00000000000..cf89c761fca --- /dev/null +++ b/docs/changes/2658.bugfix.rst @@ -0,0 +1,2 @@ +Fix ``_uncert`` calculations in ``ctapipe.reco.StereoMeanCombiner``. +Add helper functions for vectorized numpy calculations as new ``ctapipe.vectorization`` module. diff --git a/src/ctapipe/reco/stereo_combination.py b/src/ctapipe/reco/stereo_combination.py index 19c561d89b2..5fc85f91e4f 100644 --- a/src/ctapipe/reco/stereo_combination.py +++ b/src/ctapipe/reco/stereo_combination.py @@ -17,6 +17,7 @@ ReconstructedEnergyContainer, ReconstructedGeometryContainer, ) +from ..vectorization import get_subarray_index, weighted_mean_std_ufunc from .utils import add_defaults_and_meta _containers = { @@ -31,38 +32,6 @@ ] -def _grouped_add(tel_data, n_array_events, indices): - """ - Calculate the group-wise sum for each array event over the - corresponding telescope events. ``indices`` is an array - that gives the index of the subarray event for each telescope event, - returned by - ``np.unique(tel_events[["obs_id", "event_id"]], return_inverse=True)`` - """ - combined_values = np.zeros(n_array_events) - np.add.at(combined_values, indices, tel_data) - return combined_values - - -def _weighted_mean_ufunc(tel_values, weights, n_array_events, indices): - # avoid numerical problems by very large or small weights - weights = weights / weights.max() - sum_prediction = _grouped_add( - tel_values * weights, - n_array_events, - indices, - ) - sum_of_weights = _grouped_add( - weights, - n_array_events, - indices, - ) - mean = np.full(n_array_events, np.nan) - valid = sum_of_weights > 0 - mean[valid] = sum_prediction[valid] / sum_of_weights[valid] - return mean - - class StereoCombiner(Component): """ Base Class for algorithms combining telescope-wise predictions to common prediction. @@ -299,26 +268,27 @@ def predict_table(self, mono_predictions: Table) -> Table: prefix = f"{self.prefix}_tel" # TODO: Integrate table quality query once its done valid = mono_predictions[f"{prefix}_is_valid"] - valid_predictions = mono_predictions[valid] - array_events, indices, multiplicity = np.unique( - mono_predictions[["obs_id", "event_id"]], - return_inverse=True, - return_counts=True, + obs_ids, event_ids, multiplicity, tel_to_array_indices = get_subarray_index( + mono_predictions ) - stereo_table = Table(array_events) + n_array_events = len(obs_ids) + stereo_table = Table({"obs_id": obs_ids, "event_id": event_ids}) # copy metadata for colname in ("obs_id", "event_id"): stereo_table[colname].description = mono_predictions[colname].description - n_array_events = len(array_events) - weights = self._calculate_weights(valid_predictions) + weights = self._calculate_weights(mono_predictions[valid]) if self.property is ReconstructionProperty.PARTICLE_TYPE: - if len(valid_predictions) > 0: - mono_predictions = valid_predictions[f"{prefix}_prediction"] - stereo_predictions = _weighted_mean_ufunc( - mono_predictions, weights, n_array_events, indices[valid] + if np.sum(valid) > 0: + stereo_predictions, _ = weighted_mean_std_ufunc( + mono_predictions[f"{prefix}_prediction"], + valid, + n_array_events, + tel_to_array_indices, + multiplicity, + weights=weights, ) else: stereo_predictions = np.full(n_array_events, np.nan) @@ -328,29 +298,21 @@ def predict_table(self, mono_predictions: Table) -> Table: stereo_table[f"{self.prefix}_goodness_of_fit"] = np.nan elif self.property is ReconstructionProperty.ENERGY: - if len(valid_predictions) > 0: - mono_energies = valid_predictions[f"{prefix}_energy"].quantity.to_value( + if np.sum(valid) > 0: + mono_energies = mono_predictions[f"{prefix}_energy"].quantity.to_value( u.TeV ) - if self.log_target: mono_energies = np.log(mono_energies) - stereo_energy = _weighted_mean_ufunc( + stereo_energy, std = weighted_mean_std_ufunc( mono_energies, - weights, + valid, n_array_events, - indices[valid], + tel_to_array_indices, + multiplicity, + weights=weights, ) - variance = _weighted_mean_ufunc( - (mono_energies - np.repeat(stereo_energy, multiplicity)[valid]) - ** 2, - weights, - n_array_events, - indices[valid], - ) - std = np.sqrt(variance) - if self.log_target: stereo_energy = np.exp(stereo_energy) std = np.exp(std) @@ -369,22 +331,37 @@ def predict_table(self, mono_predictions: Table) -> Table: stereo_table[f"{self.prefix}_goodness_of_fit"] = np.nan elif self.property is ReconstructionProperty.GEOMETRY: - if len(valid_predictions) > 0: + if np.sum(valid) > 0: coord = AltAz( - alt=valid_predictions[f"{prefix}_alt"], - az=valid_predictions[f"{prefix}_az"], + alt=mono_predictions[f"{prefix}_alt"], + az=mono_predictions[f"{prefix}_az"], ) # https://en.wikipedia.org/wiki/Von_Mises%E2%80%93Fisher_distribution#Mean_direction mono_x, mono_y, mono_z = coord.cartesian.get_xyz() - stereo_x = _weighted_mean_ufunc( - mono_x, weights, n_array_events, indices[valid] + stereo_x, _ = weighted_mean_std_ufunc( + mono_x, + valid, + n_array_events, + tel_to_array_indices, + multiplicity, + weights=weights, ) - stereo_y = _weighted_mean_ufunc( - mono_y, weights, n_array_events, indices[valid] + stereo_y, _ = weighted_mean_std_ufunc( + mono_y, + valid, + n_array_events, + tel_to_array_indices, + multiplicity, + weights=weights, ) - stereo_z = _weighted_mean_ufunc( - mono_z, weights, n_array_events, indices[valid] + stereo_z, _ = weighted_mean_std_ufunc( + mono_z, + valid, + n_array_events, + tel_to_array_indices, + multiplicity, + weights=weights, ) mean_cartesian = CartesianRepresentation( @@ -425,7 +402,9 @@ def predict_table(self, mono_predictions: Table) -> Table: tel_ids = [[] for _ in range(n_array_events)] - for index, tel_id in zip(indices[valid], valid_predictions["tel_id"]): + for index, tel_id in zip( + tel_to_array_indices[valid], mono_predictions["tel_id"][valid] + ): tel_ids[index].append(tel_id) stereo_table[f"{self.prefix}_telescopes"] = tel_ids diff --git a/src/ctapipe/reco/tests/test_stereo_combination.py b/src/ctapipe/reco/tests/test_stereo_combination.py index 7187262af27..63ced2a980c 100644 --- a/src/ctapipe/reco/tests/test_stereo_combination.py +++ b/src/ctapipe/reco/tests/test_stereo_combination.py @@ -84,8 +84,18 @@ def test_predict_mean_energy(weights, mono_table): assert_array_equal(stereo["event_id"], np.array([1, 2, 1])) if weights == "intensity": assert_array_equal(stereo["dummy_energy"], [7, 0.5, np.nan] * u.TeV) + assert_allclose( + stereo["dummy_energy_uncert"].quantity, + [4.242641, 0, np.nan] * u.TeV, + atol=1e-7, + ) elif weights == "none": assert_array_equal(stereo["dummy_energy"], [5, 0.5, np.nan] * u.TeV) + assert_allclose( + stereo["dummy_energy_uncert"].quantity, + [3.741657, 0, np.nan] * u.TeV, + atol=1e-7, + ) assert_array_equal(stereo["dummy_telescopes"][0], np.array([1, 2, 3])) assert_array_equal(stereo["dummy_telescopes"][1], 5) diff --git a/src/ctapipe/tools/apply_models.py b/src/ctapipe/tools/apply_models.py index 6ec237006e0..3caff55fed1 100644 --- a/src/ctapipe/tools/apply_models.py +++ b/src/ctapipe/tools/apply_models.py @@ -1,6 +1,7 @@ """ Tool to apply machine learning models in bulk (as opposed to event by event). """ + import numpy as np import tables from astropy.table import Table, join, vstack @@ -9,9 +10,8 @@ from ctapipe.core.tool import Tool from ctapipe.core.traits import Bool, Integer, List, Path, classes_with_traits, flag from ctapipe.io import HDF5Merger, TableLoader, write_table -from ctapipe.io.astropy_helpers import read_table +from ctapipe.io.astropy_helpers import join_allow_empty, read_table from ctapipe.io.tableio import TelListToMaskTransform -from ctapipe.io.tableloader import _join_subarray_events from ctapipe.reco import Reconstructor __all__ = [ @@ -218,6 +218,23 @@ def _apply(self, reconstructor, tel_tables, start, stop): def _combine(self, reconstructor, tel_tables, start, stop): stereo_table = vstack(list(tel_tables.values())) + # stacking the single telescope tables and joining + # potentially changes the order of the subarray events. + # to ensure events are stored in the correct order, + # we resort to trigger table order + trigger = read_table( + self.h5file, "/dl1/event/subarray/trigger", start=start, stop=stop + )[["obs_id", "event_id"]] + trigger["__sort_index__"] = np.arange(len(trigger)) + + stereo_table = join_allow_empty( + stereo_table, + trigger, + keys=["obs_id", "event_id"], + join_type="left", + ) + stereo_table.sort("__sort_index__") + combiner = reconstructor.stereo_combiner stereo_predictions = combiner.predict_table(stereo_table) del stereo_table @@ -230,19 +247,6 @@ def _combine(self, reconstructor, tel_tables, start, stop): stereo_predictions[c.name] = np.array([trafo(r) for r in c]) stereo_predictions[c.name].description = c.description - # stacking the single telescope tables and joining - # potentially changes the order of the subarray events. - # to ensure events are stored in the correct order, - # we resort to trigger table order - trigger = read_table( - self.h5file, "/dl1/event/subarray/trigger", start=start, stop=stop - )[["obs_id", "event_id"]] - trigger["__sort_index__"] = np.arange(len(trigger)) - - stereo_predictions = _join_subarray_events(trigger, stereo_predictions) - stereo_predictions.sort("__sort_index__") - del stereo_predictions["__sort_index__"] - write_table( stereo_predictions, self.output_path, diff --git a/src/ctapipe/vectorization/__init__.py b/src/ctapipe/vectorization/__init__.py new file mode 100644 index 00000000000..b31ac7d1c6a --- /dev/null +++ b/src/ctapipe/vectorization/__init__.py @@ -0,0 +1,6 @@ +from .aggregate import get_subarray_index, weighted_mean_std_ufunc + +__all__ = [ + "get_subarray_index", + "weighted_mean_std_ufunc", +] diff --git a/src/ctapipe/vectorization/aggregate.py b/src/ctapipe/vectorization/aggregate.py new file mode 100644 index 00000000000..2c5ffc3d3b1 --- /dev/null +++ b/src/ctapipe/vectorization/aggregate.py @@ -0,0 +1,145 @@ +"""Helper functions for vectorizing numpy operations.""" + +import numpy as np +from numba import njit, uint64 + +__all__ = ["get_subarray_index", "weighted_mean_std_ufunc"] + + +@njit +def _get_subarray_index(obs_ids, event_ids): + n_tel_events = len(obs_ids) + idx = np.zeros(n_tel_events, dtype=uint64) + current_idx = 0 + subarray_obs_index = [] + subarray_event_index = [] + multiplicities = [] + multiplicity = 0 + + if n_tel_events > 0: + subarray_obs_index.append(obs_ids[0]) + subarray_event_index.append(event_ids[0]) + multiplicity += 1 + + for i in range(1, n_tel_events): + if obs_ids[i] != obs_ids[i - 1] or event_ids[i] != event_ids[i - 1]: + # append to subarray events + multiplicities.append(multiplicity) + subarray_obs_index.append(obs_ids[i]) + subarray_event_index.append(event_ids[i]) + # reset state + current_idx += 1 + multiplicity = 0 + + multiplicity += 1 + idx[i] = current_idx + + # Append multiplicity of the last subarray event + if n_tel_events > 0: + multiplicities.append(multiplicity) + + return ( + np.asarray(subarray_obs_index), + np.asarray(subarray_event_index), + np.asarray(multiplicities), + idx, + ) + + +def get_subarray_index(tel_table): + """ + Get the obs_ids and event_ids of all subarray events contained + in a table of telescope events, their multiplicity and an array + giving the index of the subarray event for each telescope event. + This requires the telescope events to be SORTED by their corresponding + subarray events (meaning by ``["obs_id", "event_id"]``). + + Parameters + ---------- + tel_table: astropy.table.Table + table with telescope events as rows + + Returns + ------- + Tuple(np.ndarray, np.ndarray, np.ndarray, np.ndarray) + obs_ids of subarray events, event_ids of subarray events, + multiplicity of subarray events, index of the subarray event + for each telescope event + """ + obs_idx = tel_table["obs_id"] + event_idx = tel_table["event_id"] + return _get_subarray_index(obs_idx, event_idx) + + +def _grouped_add(tel_data, n_array_events, indices): + """ + Calculate the group-wise sum for each array event over the + corresponding telescope events. ``indices`` is an array + that gives the index of the subarray event for each telescope event. + """ + combined_values = np.zeros(n_array_events) + np.add.at(combined_values, indices, tel_data) + return combined_values + + +def weighted_mean_std_ufunc( + tel_values, + valid_tel, + n_array_events, + indices, + multiplicity, + weights=np.array([1]), +): + """ + Calculate the weighted mean and standard deviation for each array event + over the corresponding telescope events. + + Parameters + ---------- + tel_values: np.ndarray + values for each telescope event + valid_tel: array-like + boolean mask giving the valid values of ``tel_values`` + n_array_events: int + number of array events with corresponding telescope events in ``tel_values`` + indices: np.ndarray + index of the subarray event for each telescope event, returned as + the fourth return value of ``get_subarray_index`` + multiplicity: np.ndarray + multiplicity of the subarray events in the same order as the order of + subarray events in ``indices`` + weights: np.ndarray + weights used for averaging (equal/no weights are used by default) + + Returns + ------- + Tuple(np.ndarray, np.ndarray) + weighted mean and standard deviation for each array event + """ + # avoid numerical problems by very large or small weights + weights = weights / weights.max() + tel_values = tel_values[valid_tel] + indices = indices[valid_tel] + + sum_prediction = _grouped_add( + tel_values * weights, + n_array_events, + indices, + ) + sum_of_weights = _grouped_add( + weights, + n_array_events, + indices, + ) + mean = np.full(n_array_events, np.nan) + valid = sum_of_weights > 0 + mean[valid] = sum_prediction[valid] / sum_of_weights[valid] + + sum_sq_residulas = _grouped_add( + (tel_values - np.repeat(mean, multiplicity)[valid_tel]) ** 2 * weights, + n_array_events, + indices, + ) + variance = np.full(n_array_events, np.nan) + variance[valid] = sum_sq_residulas[valid] / sum_of_weights[valid] + return mean, np.sqrt(variance) diff --git a/src/ctapipe/vectorization/tests/test_aggregate.py b/src/ctapipe/vectorization/tests/test_aggregate.py new file mode 100644 index 00000000000..14c3f578c4f --- /dev/null +++ b/src/ctapipe/vectorization/tests/test_aggregate.py @@ -0,0 +1,50 @@ +import numpy as np +from astropy.table import Table + +from ctapipe.io import TableLoader, read_table +from ctapipe.io.tests.test_table_loader import check_equal_array_event_order + + +def test_get_subarray_index(dl1_parameters_file): + from ctapipe.vectorization import get_subarray_index + + opts = dict(simulated=False, true_parameters=False, dl2=False, pointing=False) + with TableLoader(dl1_parameters_file, **opts) as loader: + tel_events = loader.read_telescope_events() + + subarray_obs_ids, subarray_event_ids, _, _ = get_subarray_index(tel_events) + trigger = read_table(dl1_parameters_file, "/dl1/event/subarray/trigger") + + assert len(subarray_obs_ids) == len(subarray_event_ids) + assert len(subarray_obs_ids) == len(trigger) + check_equal_array_event_order( + Table({"obs_id": subarray_obs_ids, "event_id": subarray_event_ids}), trigger + ) + + +def test_mean_std_ufunc(dl1_parameters_file): + from ctapipe.vectorization import get_subarray_index, weighted_mean_std_ufunc + + opts = dict(simulated=False, true_parameters=False, dl2=False, pointing=False) + with TableLoader(dl1_parameters_file, **opts) as loader: + tel_events = loader.read_telescope_events() + + valid = np.isfinite(tel_events["hillas_length"]) + obs_ids, event_ids, m, tel_to_subarray_idx = get_subarray_index(tel_events) + n_array_events = len(obs_ids) + + # test only default uniform weights, + # other weights are tested in test_stereo_combination + mean, std = weighted_mean_std_ufunc( + tel_events["hillas_length"], valid, n_array_events, tel_to_subarray_idx, m + ) + + # check if result is identical with np.mean/np.std + true_mean = np.full(n_array_events, np.nan) + true_std = np.full(n_array_events, np.nan) + grouped = tel_events.group_by(["obs_id", "event_id"]) + true_mean = grouped["hillas_length"].groups.aggregate(np.nanmean) + true_std = grouped["hillas_length"].groups.aggregate(np.nanstd) + + assert np.allclose(mean, true_mean, equal_nan=True) + assert np.allclose(std, true_std, equal_nan=True) From 2e83a98779b2496ebdd75474ff6ec48a68e2123c Mon Sep 17 00:00:00 2001 From: Lukas Beiske Date: Tue, 26 Nov 2024 12:32:25 +0100 Subject: [PATCH 2/8] Move helper functions into ctapipe.reco submodule --- docs/api-reference/vectorization/index.rst | 20 ------------------- src/ctapipe/reco/stereo_combination.py | 2 +- .../telecope_event_handling.py} | 2 +- .../tests/test_telescope_event_handling.py} | 7 +++++-- src/ctapipe/vectorization/__init__.py | 6 ------ 5 files changed, 7 insertions(+), 30 deletions(-) delete mode 100644 docs/api-reference/vectorization/index.rst rename src/ctapipe/{vectorization/aggregate.py => reco/telecope_event_handling.py} (98%) rename src/ctapipe/{vectorization/tests/test_aggregate.py => reco/tests/test_telescope_event_handling.py} (90%) delete mode 100644 src/ctapipe/vectorization/__init__.py diff --git a/docs/api-reference/vectorization/index.rst b/docs/api-reference/vectorization/index.rst deleted file mode 100644 index 640a28b2eb0..00000000000 --- a/docs/api-reference/vectorization/index.rst +++ /dev/null @@ -1,20 +0,0 @@ -.. _vectorization: - -******************************************************************************* -Helper Functions for Vectorized Operations on Tables (`~ctapipe.vectorization`) -******************************************************************************* - -.. currentmodule:: ctapipe.vectorization - -This module contains helper functions for performing vectorized operations -on tables, as used in e.g. :ref:`reco_stereo_combination`. - - -Reference/API -============= - -.. automodapi:: ctapipe.vectorization - :no-inheritance-diagram: - -.. automodapi:: ctapipe.vectorization.aggregate - :no-inheritance-diagram: diff --git a/src/ctapipe/reco/stereo_combination.py b/src/ctapipe/reco/stereo_combination.py index 5fc85f91e4f..f35b2c7d794 100644 --- a/src/ctapipe/reco/stereo_combination.py +++ b/src/ctapipe/reco/stereo_combination.py @@ -17,7 +17,7 @@ ReconstructedEnergyContainer, ReconstructedGeometryContainer, ) -from ..vectorization import get_subarray_index, weighted_mean_std_ufunc +from .telecope_event_handling import get_subarray_index, weighted_mean_std_ufunc from .utils import add_defaults_and_meta _containers = { diff --git a/src/ctapipe/vectorization/aggregate.py b/src/ctapipe/reco/telecope_event_handling.py similarity index 98% rename from src/ctapipe/vectorization/aggregate.py rename to src/ctapipe/reco/telecope_event_handling.py index 2c5ffc3d3b1..2fc0e3233cd 100644 --- a/src/ctapipe/vectorization/aggregate.py +++ b/src/ctapipe/reco/telecope_event_handling.py @@ -1,4 +1,4 @@ -"""Helper functions for vectorizing numpy operations.""" +"""Helper functions for array-event-wise aggregation of telescope events.""" import numpy as np from numba import njit, uint64 diff --git a/src/ctapipe/vectorization/tests/test_aggregate.py b/src/ctapipe/reco/tests/test_telescope_event_handling.py similarity index 90% rename from src/ctapipe/vectorization/tests/test_aggregate.py rename to src/ctapipe/reco/tests/test_telescope_event_handling.py index 14c3f578c4f..c238029fd90 100644 --- a/src/ctapipe/vectorization/tests/test_aggregate.py +++ b/src/ctapipe/reco/tests/test_telescope_event_handling.py @@ -6,7 +6,7 @@ def test_get_subarray_index(dl1_parameters_file): - from ctapipe.vectorization import get_subarray_index + from ctapipe.reco.telecope_event_handling import get_subarray_index opts = dict(simulated=False, true_parameters=False, dl2=False, pointing=False) with TableLoader(dl1_parameters_file, **opts) as loader: @@ -23,7 +23,10 @@ def test_get_subarray_index(dl1_parameters_file): def test_mean_std_ufunc(dl1_parameters_file): - from ctapipe.vectorization import get_subarray_index, weighted_mean_std_ufunc + from ctapipe.reco.telecope_event_handling import ( + get_subarray_index, + weighted_mean_std_ufunc, + ) opts = dict(simulated=False, true_parameters=False, dl2=False, pointing=False) with TableLoader(dl1_parameters_file, **opts) as loader: diff --git a/src/ctapipe/vectorization/__init__.py b/src/ctapipe/vectorization/__init__.py deleted file mode 100644 index b31ac7d1c6a..00000000000 --- a/src/ctapipe/vectorization/__init__.py +++ /dev/null @@ -1,6 +0,0 @@ -from .aggregate import get_subarray_index, weighted_mean_std_ufunc - -__all__ = [ - "get_subarray_index", - "weighted_mean_std_ufunc", -] From ca1bcfdc2e22ff3cd7360dfaaca8c30bb16fa96a Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Mon, 2 Dec 2024 14:09:22 +0100 Subject: [PATCH 3/8] sum -> count_nonzero --- src/ctapipe/reco/stereo_combination.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/ctapipe/reco/stereo_combination.py b/src/ctapipe/reco/stereo_combination.py index f35b2c7d794..a90c26069ff 100644 --- a/src/ctapipe/reco/stereo_combination.py +++ b/src/ctapipe/reco/stereo_combination.py @@ -281,7 +281,7 @@ def predict_table(self, mono_predictions: Table) -> Table: weights = self._calculate_weights(mono_predictions[valid]) if self.property is ReconstructionProperty.PARTICLE_TYPE: - if np.sum(valid) > 0: + if np.count_nonzero(valid) > 0: stereo_predictions, _ = weighted_mean_std_ufunc( mono_predictions[f"{prefix}_prediction"], valid, @@ -298,7 +298,7 @@ def predict_table(self, mono_predictions: Table) -> Table: stereo_table[f"{self.prefix}_goodness_of_fit"] = np.nan elif self.property is ReconstructionProperty.ENERGY: - if np.sum(valid) > 0: + if np.count_nonzero(valid) > 0: mono_energies = mono_predictions[f"{prefix}_energy"].quantity.to_value( u.TeV ) @@ -331,7 +331,7 @@ def predict_table(self, mono_predictions: Table) -> Table: stereo_table[f"{self.prefix}_goodness_of_fit"] = np.nan elif self.property is ReconstructionProperty.GEOMETRY: - if np.sum(valid) > 0: + if np.count_nonzero(valid) > 0: coord = AltAz( alt=mono_predictions[f"{prefix}_alt"], az=mono_predictions[f"{prefix}_az"], From 8b4e29baf8a09708a172240e1080c759dca24c90 Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Mon, 2 Dec 2024 14:10:56 +0100 Subject: [PATCH 4/8] Fix typo in module name --- src/ctapipe/coordinates/tests/test_impact_distance.py | 4 ++-- src/ctapipe/reco/stereo_combination.py | 2 +- ...telecope_event_handling.py => telescope_event_handling.py} | 0 src/ctapipe/reco/tests/test_telescope_event_handling.py | 4 ++-- 4 files changed, 5 insertions(+), 5 deletions(-) rename src/ctapipe/reco/{telecope_event_handling.py => telescope_event_handling.py} (100%) diff --git a/src/ctapipe/coordinates/tests/test_impact_distance.py b/src/ctapipe/coordinates/tests/test_impact_distance.py index c5093361cee..fb7b044c15c 100644 --- a/src/ctapipe/coordinates/tests/test_impact_distance.py +++ b/src/ctapipe/coordinates/tests/test_impact_distance.py @@ -75,7 +75,7 @@ def test_shower_impact_distance(reference_location): assert np.allclose(impact_distances, [0, 1, 2] * u.m) # alt=0 az=0 should be aligned to x-axis (north in ground frame) - # therefore the distances should also be just the y-offset of the telecope + # therefore the distances should also be just the y-offset of the telescope shower_geom = ReconstructedGeometryContainer( core_x=0 * u.m, core_y=0 * u.m, alt=0 * u.deg, az=0 * u.deg ) @@ -83,7 +83,7 @@ def test_shower_impact_distance(reference_location): assert np.allclose(impact_distances, [0, 1, 2] * u.m) # alt=0 az=90 should be aligned to y-axis (east in ground frame) - # therefore the distances should also be just the x-offset of the telecope + # therefore the distances should also be just the x-offset of the telescope shower_geom = ReconstructedGeometryContainer( core_x=0 * u.m, core_y=0 * u.m, alt=0 * u.deg, az=90 * u.deg ) diff --git a/src/ctapipe/reco/stereo_combination.py b/src/ctapipe/reco/stereo_combination.py index a90c26069ff..480ed5c6ed7 100644 --- a/src/ctapipe/reco/stereo_combination.py +++ b/src/ctapipe/reco/stereo_combination.py @@ -17,7 +17,7 @@ ReconstructedEnergyContainer, ReconstructedGeometryContainer, ) -from .telecope_event_handling import get_subarray_index, weighted_mean_std_ufunc +from .telescope_event_handling import get_subarray_index, weighted_mean_std_ufunc from .utils import add_defaults_and_meta _containers = { diff --git a/src/ctapipe/reco/telecope_event_handling.py b/src/ctapipe/reco/telescope_event_handling.py similarity index 100% rename from src/ctapipe/reco/telecope_event_handling.py rename to src/ctapipe/reco/telescope_event_handling.py diff --git a/src/ctapipe/reco/tests/test_telescope_event_handling.py b/src/ctapipe/reco/tests/test_telescope_event_handling.py index c238029fd90..86f0a112345 100644 --- a/src/ctapipe/reco/tests/test_telescope_event_handling.py +++ b/src/ctapipe/reco/tests/test_telescope_event_handling.py @@ -6,7 +6,7 @@ def test_get_subarray_index(dl1_parameters_file): - from ctapipe.reco.telecope_event_handling import get_subarray_index + from ctapipe.reco.telescope_event_handling import get_subarray_index opts = dict(simulated=False, true_parameters=False, dl2=False, pointing=False) with TableLoader(dl1_parameters_file, **opts) as loader: @@ -23,7 +23,7 @@ def test_get_subarray_index(dl1_parameters_file): def test_mean_std_ufunc(dl1_parameters_file): - from ctapipe.reco.telecope_event_handling import ( + from ctapipe.reco.telescope_event_handling import ( get_subarray_index, weighted_mean_std_ufunc, ) From c0ba9862845e8643c489af7ceef21c138c586f0e Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Mon, 2 Dec 2024 14:23:46 +0100 Subject: [PATCH 5/8] Better docstring --- src/ctapipe/reco/telescope_event_handling.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/ctapipe/reco/telescope_event_handling.py b/src/ctapipe/reco/telescope_event_handling.py index 2fc0e3233cd..c2e7a5c4d0d 100644 --- a/src/ctapipe/reco/telescope_event_handling.py +++ b/src/ctapipe/reco/telescope_event_handling.py @@ -48,11 +48,14 @@ def _get_subarray_index(obs_ids, event_ids): def get_subarray_index(tel_table): """ - Get the obs_ids and event_ids of all subarray events contained + Get the subarray-event-wise information from a table of telescope events. + + Extract obs_ids and event_ids of all subarray events contained in a table of telescope events, their multiplicity and an array giving the index of the subarray event for each telescope event. - This requires the telescope events to be SORTED by their corresponding - subarray events (meaning by ``["obs_id", "event_id"]``). + + This requires the telescope events of one subarray event to be come + in a single block. Parameters ---------- From 3e7269175815e2f83533186a4085c563a953554f Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Mon, 2 Dec 2024 14:29:34 +0100 Subject: [PATCH 6/8] Remove unneeded argument from weighted_mean_std --- src/ctapipe/reco/stereo_combination.py | 5 ----- src/ctapipe/reco/telescope_event_handling.py | 4 +--- src/ctapipe/reco/tests/test_telescope_event_handling.py | 7 ++++--- 3 files changed, 5 insertions(+), 11 deletions(-) diff --git a/src/ctapipe/reco/stereo_combination.py b/src/ctapipe/reco/stereo_combination.py index 480ed5c6ed7..6261af07476 100644 --- a/src/ctapipe/reco/stereo_combination.py +++ b/src/ctapipe/reco/stereo_combination.py @@ -285,7 +285,6 @@ def predict_table(self, mono_predictions: Table) -> Table: stereo_predictions, _ = weighted_mean_std_ufunc( mono_predictions[f"{prefix}_prediction"], valid, - n_array_events, tel_to_array_indices, multiplicity, weights=weights, @@ -308,7 +307,6 @@ def predict_table(self, mono_predictions: Table) -> Table: stereo_energy, std = weighted_mean_std_ufunc( mono_energies, valid, - n_array_events, tel_to_array_indices, multiplicity, weights=weights, @@ -342,7 +340,6 @@ def predict_table(self, mono_predictions: Table) -> Table: stereo_x, _ = weighted_mean_std_ufunc( mono_x, valid, - n_array_events, tel_to_array_indices, multiplicity, weights=weights, @@ -350,7 +347,6 @@ def predict_table(self, mono_predictions: Table) -> Table: stereo_y, _ = weighted_mean_std_ufunc( mono_y, valid, - n_array_events, tel_to_array_indices, multiplicity, weights=weights, @@ -358,7 +354,6 @@ def predict_table(self, mono_predictions: Table) -> Table: stereo_z, _ = weighted_mean_std_ufunc( mono_z, valid, - n_array_events, tel_to_array_indices, multiplicity, weights=weights, diff --git a/src/ctapipe/reco/telescope_event_handling.py b/src/ctapipe/reco/telescope_event_handling.py index c2e7a5c4d0d..5e9da9621e8 100644 --- a/src/ctapipe/reco/telescope_event_handling.py +++ b/src/ctapipe/reco/telescope_event_handling.py @@ -88,7 +88,6 @@ def _grouped_add(tel_data, n_array_events, indices): def weighted_mean_std_ufunc( tel_values, valid_tel, - n_array_events, indices, multiplicity, weights=np.array([1]), @@ -103,8 +102,6 @@ def weighted_mean_std_ufunc( values for each telescope event valid_tel: array-like boolean mask giving the valid values of ``tel_values`` - n_array_events: int - number of array events with corresponding telescope events in ``tel_values`` indices: np.ndarray index of the subarray event for each telescope event, returned as the fourth return value of ``get_subarray_index`` @@ -119,6 +116,7 @@ def weighted_mean_std_ufunc( Tuple(np.ndarray, np.ndarray) weighted mean and standard deviation for each array event """ + n_array_events = len(multiplicity) # avoid numerical problems by very large or small weights weights = weights / weights.max() tel_values = tel_values[valid_tel] diff --git a/src/ctapipe/reco/tests/test_telescope_event_handling.py b/src/ctapipe/reco/tests/test_telescope_event_handling.py index 86f0a112345..6e74fcdcc08 100644 --- a/src/ctapipe/reco/tests/test_telescope_event_handling.py +++ b/src/ctapipe/reco/tests/test_telescope_event_handling.py @@ -33,13 +33,14 @@ def test_mean_std_ufunc(dl1_parameters_file): tel_events = loader.read_telescope_events() valid = np.isfinite(tel_events["hillas_length"]) - obs_ids, event_ids, m, tel_to_subarray_idx = get_subarray_index(tel_events) - n_array_events = len(obs_ids) + + _, _, multiplicity, tel_to_subarray_idx = get_subarray_index(tel_events) + n_array_events = len(multiplicity) # test only default uniform weights, # other weights are tested in test_stereo_combination mean, std = weighted_mean_std_ufunc( - tel_events["hillas_length"], valid, n_array_events, tel_to_subarray_idx, m + tel_events["hillas_length"], valid, tel_to_subarray_idx, multiplicity ) # check if result is identical with np.mean/np.std From c380912dbecd238b19514c8e237b4125e3594b29 Mon Sep 17 00:00:00 2001 From: Lukas Beiske Date: Mon, 2 Dec 2024 14:37:02 +0100 Subject: [PATCH 7/8] Remove unused variables in test --- src/ctapipe/reco/tests/test_telescope_event_handling.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/ctapipe/reco/tests/test_telescope_event_handling.py b/src/ctapipe/reco/tests/test_telescope_event_handling.py index 6e74fcdcc08..1faa0f4f65c 100644 --- a/src/ctapipe/reco/tests/test_telescope_event_handling.py +++ b/src/ctapipe/reco/tests/test_telescope_event_handling.py @@ -35,7 +35,6 @@ def test_mean_std_ufunc(dl1_parameters_file): valid = np.isfinite(tel_events["hillas_length"]) _, _, multiplicity, tel_to_subarray_idx = get_subarray_index(tel_events) - n_array_events = len(multiplicity) # test only default uniform weights, # other weights are tested in test_stereo_combination @@ -44,8 +43,6 @@ def test_mean_std_ufunc(dl1_parameters_file): ) # check if result is identical with np.mean/np.std - true_mean = np.full(n_array_events, np.nan) - true_std = np.full(n_array_events, np.nan) grouped = tel_events.group_by(["obs_id", "event_id"]) true_mean = grouped["hillas_length"].groups.aggregate(np.nanmean) true_std = grouped["hillas_length"].groups.aggregate(np.nanstd) From 9c74e4d76d6e81a1de1c74ec56a14fb3daabfc20 Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Mon, 2 Dec 2024 17:11:11 +0100 Subject: [PATCH 8/8] Fix changelog entry --- docs/changes/2658.bugfix.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/changes/2658.bugfix.rst b/docs/changes/2658.bugfix.rst index cf89c761fca..40b1e8239e6 100644 --- a/docs/changes/2658.bugfix.rst +++ b/docs/changes/2658.bugfix.rst @@ -1,2 +1,2 @@ Fix ``_uncert`` calculations in ``ctapipe.reco.StereoMeanCombiner``. -Add helper functions for vectorized numpy calculations as new ``ctapipe.vectorization`` module. +Add helper functions for vectorized numpy calculations as new ``ctapipe.reco.telescope_event_handling`` module.