diff --git a/docs/changes/2731.feature.rst b/docs/changes/2731.feature.rst new file mode 100644 index 00000000000..ffca00e7784 --- /dev/null +++ b/docs/changes/2731.feature.rst @@ -0,0 +1,11 @@ +Added a new ``StereoDispCombiner`` implementing DISP-based stereo direction +reconstruction. +The algorithm is adapted and generalized from the implementations in +*magic-cta-pipe* and *EventDisplay* and resolves the head–tail ambiguity by +evaluating all DISP sign combinations and combining multiple telescope +combinations via weighted means. + +The reconstruction supports optional pre-selection of telescopes using +``n_best_tels`` and rejection of nearly parallel multiplicity-2 events via +``min_ang_diff``. +Only geometry reconstruction is supported. diff --git a/src/ctapipe/image/statistics.py b/src/ctapipe/image/statistics.py index f0d94a81711..332fe00fbac 100644 --- a/src/ctapipe/image/statistics.py +++ b/src/ctapipe/image/statistics.py @@ -104,7 +104,7 @@ def descriptive_statistics( ) -@njit +@njit(cache=not CTAPIPE_DISABLE_NUMBA_CACHE) def n_largest(n, array): """return the n largest values of an array""" return nlargest(n, array) diff --git a/src/ctapipe/reco/__init__.py b/src/ctapipe/reco/__init__.py index 267b68b4b5c..464b1280426 100644 --- a/src/ctapipe/reco/__init__.py +++ b/src/ctapipe/reco/__init__.py @@ -17,7 +17,7 @@ EnergyRegressor, ParticleClassifier, ) -from .stereo_combination import StereoCombiner, StereoMeanCombiner +from .stereo_combination import StereoCombiner, StereoDispCombiner, StereoMeanCombiner __all__ = [ "Reconstructor", @@ -32,5 +32,6 @@ "DispReconstructor", "StereoCombiner", "StereoMeanCombiner", + "StereoDispCombiner", "CrossValidator", ] diff --git a/src/ctapipe/reco/stereo_combination.py b/src/ctapipe/reco/stereo_combination.py index b28c93713ca..1c1b5929ce0 100644 --- a/src/ctapipe/reco/stereo_combination.py +++ b/src/ctapipe/reco/stereo_combination.py @@ -6,8 +6,17 @@ from astropy.table import Table from traitlets import UseEnum -from ctapipe.core import Component, Container -from ctapipe.core.traits import Bool, CaselessStrEnum, Unicode +from ctapipe.containers import ImageParametersContainer +from ctapipe.core import Component +from ctapipe.core.traits import ( + Bool, + CaselessStrEnum, + Float, + Integer, + TraitError, + Unicode, +) +from ctapipe.image.statistics import arg_n_largest from ctapipe.reco.reconstructor import ReconstructionProperty from ..compat import COPY_IF_NEEDED @@ -17,7 +26,20 @@ ReconstructedEnergyContainer, ReconstructedGeometryContainer, ) -from .telescope_event_handling import get_subarray_index, weighted_mean_std_ufunc +from .preprocessing import telescope_to_horizontal +from .reconstructor import StereoQualityQuery +from .telescope_event_handling import ( + calc_combs_min_distances, + calc_fov_lon_lat, + check_ang_diff, + create_combs_array, + fill_lower_multiplicities, + get_combinations, + get_index_combs, + get_subarray_index, + valid_tels_of_multi, + weighted_mean_std_ufunc, +) from .utils import add_defaults_and_meta _containers = { @@ -29,71 +51,73 @@ __all__ = [ "StereoCombiner", "StereoMeanCombiner", + "StereoDispCombiner", ] class StereoCombiner(Component): """ - Base Class for algorithms combining telescope-wise predictions to common prediction. - """ + Base class for algorithms combining telescope-wise predictions + to common predictions. - prefix = Unicode( - default_value="", - help="Prefix to be added to the output container / column names.", - ).tag(config=True) + A StereoCombiner defines the interface for transforming per-telescope + DL2 predictions (energy, direction, particle type) into a unified + stereo estimate. Subclasses implement the actual combination strategy. - property = UseEnum( - ReconstructionProperty, - help="Which property is being combined.", - ).tag(config=True) + Two usage modes are supported: - @abstractmethod - def __call__(self, event: ArrayEventContainer) -> None: - """ - Fill event container with stereo predictions. - """ + 1. **Event-wise combination** via :meth:`__call__` + Operates directly on an :class:`~ctapipe.containers.ArrayEventContainer` + containing DL1 and DL2 mono predictions for a single event, and writes + the stereo-level result into ``event.dl2.stereo``. - @abstractmethod - def predict_table(self, mono_predictions: Table) -> Table: - """ - Constructs stereo predictions from a table of - telescope events. - """ + 2. **Table-wise combination** via :meth:`predict_table` + Takes a table of DL2 mono predictions (one row per telescope-event) + and constructs an output table with one row per subarray event containing + the stereo predictions. - -class StereoMeanCombiner(StereoCombiner): - """ - Calculate array-event prediction as (weighted) mean of telescope-wise predictions. + Subclasses must implement both methods, as well as any additional logic + required for combining the chosen :class:`~ctapipe.reco.ReconstructionProperty`. """ + classes = [StereoQualityQuery] + weights = CaselessStrEnum( ["none", "intensity", "aspect-weighted-intensity"], default_value="none", help=( - "What kind of weights to use. Options: ``none``, ``intensity``, ``aspect-weighted-intensity``." + "Weighting scheme used to combine telescope-wise mono predictions:\n\n" + "- ``none``: All telescopes contribute equally (unweighted mean).\n" + "- ``intensity``: Contributions are weighted by the image intensity, " + "giving higher weight to brighter images.\n" + "- ``aspect-weighted-intensity``: Contributions are weighted by the " + "image intensity multiplied by the image aspect ratio (length/width), " + "favoring bright and well-elongated images.\n\n" + "The selected weighting is applied consistently to all supported " + "reconstruction properties (e.g. geometry, energy, particle type)." ), ).tag(config=True) - log_target = Bool( - False, - help="If true, calculate exp(mean(log(values))).", + prefix = Unicode( + default_value="", + help="Prefix to be added to the output container / column names.", + ).tag(config=True) + + property = UseEnum( + ReconstructionProperty, + help=( + "Reconstruction property to be combined (e.g. ENERGY, GEOMETRY, " + "PARTICLE_TYPE). Subclasses may support only a subset." + ), ).tag(config=True) def __init__(self, *args, **kwargs): super().__init__(*args, **kwargs) - self.supported = { - ReconstructionProperty.ENERGY, - ReconstructionProperty.GEOMETRY, - ReconstructionProperty.PARTICLE_TYPE, - } - if self.property not in self.supported: - raise NotImplementedError( - f"Combination of {self.property} not implemented in {self.__class__.__name__}" - ) + self.quality_query = StereoQualityQuery(parent=self) def _calculate_weights(self, data): - if isinstance(data, Container): + if isinstance(data, ImageParametersContainer): if self.weights == "intensity": return data.hillas.intensity @@ -116,9 +140,80 @@ def _calculate_weights(self, data): return np.ones(len(data)) raise TypeError( - "Dl1 data needs to be provided in the form of a container or astropy.table.Table" + "Dl1 data needs to be provided in the form of an ImageParametersContainer " + "or astropy.table.Table." ) + @abstractmethod + def __call__(self, event: ArrayEventContainer) -> None: + """ + Compute the stereo prediction for a single subarray event. + + Parameters + ---------- + event : ArrayEventContainer + Subarray event containing DL1 parameters and per-telescope DL2 predictions. + Implementations must write the combined stereo prediction into + ``event.dl2.stereo`` under the configured ``prefix``. + + Notes + ----- + - This method modifies the subarray event container *in place*. + """ + + @abstractmethod + def predict_table(self, mono_predictions: Table) -> Table: + """ + Construct stereo predictions from a table of mono DL2 predictions. + + Parameters + ---------- + mono_predictions : astropy.table.Table + Table containing one row per telescope-event and the + DL2 output corresponding to the configured reconstruction property. + + Returns + ------- + stereo_table : astropy.table.Table + A table with one row per subarray event, containing the combined + stereo predictions. + """ + + +class StereoMeanCombiner(StereoCombiner): + """ + Combine telescope-wise mono reconstructions using a (weighted) mean. + + This implementation supports combining energy, geometry, and + particle-type predictions. The weighting scheme applied to the + telescope-wise inputs is controlled via the ``weights`` configuration + trait. + + If ``log_target=True`` and ``ENERGY`` is combined, the combiner computes + the geometric mean by averaging the logarithm of the energies. + + See :class:`StereoCombiner` for a description of the general interface. + """ + + log_target = Bool( + False, + help="If true, calculate exp(mean(log(values))).", + ).tag(config=True) + + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + + self.supported = { + ReconstructionProperty.ENERGY, + ReconstructionProperty.GEOMETRY, + ReconstructionProperty.PARTICLE_TYPE, + } + if self.property not in self.supported: + raise NotImplementedError( + f"Combination of {self.property} not implemented in " + "{self.__class__.__name__}." + ) + def _combine_energy(self, event): ids = [] values = [] @@ -126,14 +221,18 @@ def _combine_energy(self, event): for tel_id, dl2 in event.dl2.tel.items(): mono = dl2.energy[self.prefix] - if mono.is_valid: - values.append(mono.energy.to_value(u.TeV)) - if tel_id not in event.dl1.tel: - raise ValueError("No parameters for weighting available") - weights.append( - self._calculate_weights(event.dl1.tel[tel_id].parameters) - ) - ids.append(tel_id) + dl1_params = event.dl1.tel[tel_id].parameters + quality_checks = self.quality_query(parameters=dl1_params) + if not mono.is_valid: + continue + if not all(quality_checks): + continue + + values.append(mono.energy.to_value(u.TeV)) + if tel_id not in event.dl1.tel: + raise ValueError("No parameters for weighting available") + weights.append(self._calculate_weights(event.dl1.tel[tel_id].parameters)) + ids.append(tel_id) if len(values) > 0: weights = np.array(weights, dtype=np.float64) @@ -169,11 +268,17 @@ def _combine_classification(self, event): for tel_id, dl2 in event.dl2.tel.items(): mono = dl2.particle_type[self.prefix] - if mono.is_valid: - values.append(mono.prediction) - dl1 = event.dl1.tel[tel_id].parameters - weights.append(self._calculate_weights(dl1) if dl1 else 1) - ids.append(tel_id) + dl1_params = event.dl1.tel[tel_id].parameters + quality_checks = self.quality_query(parameters=dl1_params) + if not mono.is_valid: + continue + if not all(quality_checks): + continue + + values.append(mono.prediction) + dl1 = event.dl1.tel[tel_id].parameters + weights.append(self._calculate_weights(dl1) if dl1 else 1) + ids.append(tel_id) if len(values) > 0: mean = np.average(values, weights=weights) @@ -195,12 +300,18 @@ def _combine_altaz(self, event): for tel_id, dl2 in event.dl2.tel.items(): mono = dl2.geometry[self.prefix] - if mono.is_valid: - alt_values.append(mono.alt) - az_values.append(mono.az) - dl1 = event.dl1.tel[tel_id].parameters - weights.append(self._calculate_weights(dl1) if dl1 else 1) - ids.append(tel_id) + dl1_params = event.dl1.tel[tel_id].parameters + quality_checks = self.quality_query(parameters=dl1_params) + if not mono.is_valid: + continue + if not all(quality_checks): + continue + + alt_values.append(mono.alt) + az_values.append(mono.az) + dl1 = event.dl1.tel[tel_id].parameters + weights.append(self._calculate_weights(dl1) if dl1 else 1) + ids.append(tel_id) if len(alt_values) > 0: # by construction len(alt_values) == len(az_values) coord = AltAz(alt=alt_values, az=az_values) @@ -238,7 +349,8 @@ def _combine_altaz(self, event): def __call__(self, event: ArrayEventContainer) -> None: """ - Calculate the mean prediction for a single array event. + Implement :meth:`StereoCombiner.__call__` using a weighted mean + combination of the configured reconstruction property. """ properties = [ @@ -258,14 +370,17 @@ def __call__(self, event: ArrayEventContainer) -> None: def predict_table(self, mono_predictions: Table) -> Table: """ - Calculates the (array-)event-wise mean. - Telescope events, that are nan, get discarded. - This means you might end up with less events if - all telescope predictions of a shower are invalid. + Calculates the (subarray-)event-wise mean. + + Telescope events, that are nan, get discarded. This means + you might end up with less events if all telescope predictions + of a shower are invalid. + + See :meth:`StereoCombiner.predict_table` for the general + input/output conventions. """ prefix = f"{self.prefix}_tel" - # TODO: Integrate table quality query once its done valid = mono_predictions[f"{prefix}_is_valid"] obs_ids, event_ids, multiplicity, tel_to_array_indices = get_subarray_index( @@ -404,3 +519,525 @@ def predict_table(self, mono_predictions: Table) -> Table: stereo_table[f"{self.prefix}_telescopes"] = tel_ids add_defaults_and_meta(stereo_table, _containers[self.property], self.prefix) return stereo_table + + +class StereoDispCombiner(StereoCombiner): + """ + Stereo combination algorithm for DISP-based direction reconstruction. + + This combiner implements a generalized DISP stereo reconstruction based on + Algorithm 3 of :cite:p:`hofmann-1999-comparison` and the EventDisplay + approach, where each telescope yields two possible directions + (DISP sign = ±1). The algorithm resolves the head–tail ambiguity by + combining multiple telescope combinations and applying a weighted mean. + + **Algorithm overview** + + 0. **Pre-selection of valid telescopes (per subarray event)** + The set of telescopes participating in the stereo reconstruction is + optionally restricted before any combination step: + + - If ``n_best_tels`` is set, only the ``n_best_tels`` telescopes with + the highest weights are kept per subarray event; all other telescopes + are excluded from further processing. + - If ``min_ang_diff`` is set, subarray events of multiplicity = 2 that + have a nearly parallel main shower axis (angular difference below + ``min_ang_diff``) are rejected entirely. + + These selections reduce the effective telescope multiplicity used in + the subsequent steps. + + 1. For every remaining valid telescope, compute the two possible FoV + positions (longitude and latitude) from the Hillas centroid, the + reconstructed DISP distance, and the Hillas orientation angle + (``psi``), corresponding to DISP sign = ±1. + + 2. For each telescope combination of size ``n_tel_combinations``, evaluate + all possible DISP sign assignments and select the sign combination that + minimizes the Sum of Squared Errors (SSE) between the participating + telescopes. + + 3. Combine the resulting per-combination FoV positions using the selected + telescope weights and compute a weighted mean FoV direction for each + subarray event. + + 4. Convert the final FoV direction to horizontal coordinates (Alt/Az). + + **Handling of lower multiplicities** + + If an subarray event has an effective telescope multiplicity smaller than + ``n_tel_combinations`` (after applying ``n_best_tels`` and any angular + difference cuts), but at least two telescopes remain, the reconstruction + is performed using all available telescopes of that subevent. In this case, + just one combination is formed with a size equal to the event multiplicity + and the optimal DISP sign assignment is determined accordingly. Single- + telescope events are handled separately using the mono reconstruction. + + **Traits** + + - ``n_tel_combinations``: Size of each telescope combination (minimum 2). + Larger values increase computation time. For events with lower effective + multiplicity, all available telescopes are used. + - ``n_best_tels``: If set, restricts the reconstruction to the + ``n_best_tels`` highest-weight telescopes per event. + - ``min_ang_diff``: For multiplicity-2 events only, reject events with + nearly parallel main shower axes (difference below this angle). + - ``weights``: Telescope weights used in the per-combination and final + mean (``none``, ``intensity``, ``aspect-weighted-intensity``). + + Notes + ----- + - Only geometry (:class:`~ctapipe.reco.ReconstructionProperty.GEOMETRY`) + is supported. + - To reproduce the behavior of EventDisplay, set both + ``n_tel_combinations`` and ``n_best_tels`` to 5. + - Choosing large values of ``n_tel_combinations`` together with an + unrestricted ``n_best_tels`` can lead to a rapid increase in the number + of telescope combinations and thus to a significantly increased + processing time. + """ + + n_tel_combinations = Integer( + default_value=2, + min=2, + help=( + "Number of telescopes used per combination (minimum 2). Values " + "above 5 lead to significantly increased computation time." + ), + ).tag(config=True) + + n_best_tels = Integer( + default_value=None, + min=2, + allow_none=True, + help=( + "Select the best n_best_tels telescopes by weight for the combination. " + "Set to None to use all valid telescopes." + ), + ).tag(config=True) + + min_ang_diff = Float( + default_value=None, + min=0.0, + allow_none=True, + help=( + "Minimum angular separation (in degrees) between the reconstructed main " + "shower axes of two telescopes. This cut is applied only to subarray events " + "with exactly two participating telescopes (multiplicity = 2). Events for " + "which the angular difference between the two shower axes is smaller than " + "this value are rejected, as nearly parallel axes lead to problems solving " + "the head-tail ambiguity. Set to None to disable this cut." + ), + ).tag(config=True) + + property = UseEnum( + ReconstructionProperty, + default_value=ReconstructionProperty.GEOMETRY, + help=( + "Reconstruction property produced by this combiner (not configurable - " + "fixed to GEOMETRY)." + ), + ).tag(config=False) + + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + + if self.n_best_tels is not None and self.n_tel_combinations > self.n_best_tels: + raise TraitError( + "n_tel_combinations must be less than or equal to n_best_tels." + ) + self.supported = ReconstructionProperty.GEOMETRY + if self.property not in self.supported: + raise NotImplementedError( + f"Combination of {self.property} not implemented in {self.__class__.__name__}" + ) + + def __call__(self, event: ArrayEventContainer) -> None: + """ + Perform DISP-based stereo direction reconstruction for a single subarray + event. + + This is the entry point used for event-wise processing and stores the + resulting stereo geometry inside ``event.dl2.stereo`` under the + configured prefix. + """ + + ids, hillas_psis, fov_lon_values, fov_lat_values, weights = ( + self._collect_valid_tel_data(event) + ) + + multiplicity = len(ids) + alt = az = u.Quantity(np.nan, u.deg, copy=COPY_IF_NEEDED) + + stereo_fov_lon, stereo_fov_lat, valid = self._compute_stereo_fov( + multiplicity=multiplicity, + hillas_psis=hillas_psis, + fov_lon_values=fov_lon_values, + fov_lat_values=fov_lat_values, + weights=weights, + ) + + if valid: + alt, az = telescope_to_horizontal( + lon=stereo_fov_lon * u.deg, + lat=stereo_fov_lat * u.deg, + pointing_alt=event.monitoring.pointing.array_altitude, + pointing_az=event.monitoring.pointing.array_azimuth, + ) + + event.dl2.stereo.geometry[self.prefix] = ReconstructedGeometryContainer( + alt=alt, + az=az, + telescopes=ids, + is_valid=valid, + prefix=self.prefix, + ) + + def _collect_valid_tel_data(self, event: ArrayEventContainer): + ids = [] + hillas_psis = [] + fov_lon_values = [] + fov_lat_values = [] + weights = [] + + signs = np.array([-1, 1]) + + for tel_id, dl2 in event.dl2.tel.items(): + if not dl2.geometry[self.prefix].is_valid: + continue + + disp_reco = dl2.disp.get(self.prefix) + if disp_reco is None: + raise RuntimeError( + "No valid DISP reconstruction parameter found for " + f"prefix='{self.prefix}'). " + "Make sure to apply the DispReconstructor before using the " + "StereoDispCombiner or adapt the prefix accordingly." + ) + + dl1 = event.dl1.tel[tel_id].parameters + if not all(self.quality_query(parameters=dl1)): + continue + + hillas_fov_lon = dl1.hillas.fov_lon.to_value(u.deg) + hillas_fov_lat = dl1.hillas.fov_lat.to_value(u.deg) + hillas_psi = dl1.hillas.psi + disp = disp_reco.parameter.to_value(u.deg) + + fov_lons = hillas_fov_lon + signs * disp * np.cos(hillas_psi) + fov_lats = hillas_fov_lat + signs * disp * np.sin(hillas_psi) + + fov_lon_values.append(fov_lons) + fov_lat_values.append(fov_lats) + weights.append(self._calculate_weights(dl1) if dl1 else 1) + ids.append(tel_id) + hillas_psis.append(hillas_psi) + + return ids, hillas_psis, fov_lon_values, fov_lat_values, weights + + def _compute_stereo_fov( + self, + multiplicity: int, + hillas_psis, + fov_lon_values, + fov_lat_values, + weights, + ): + if multiplicity == 0: + return np.nan, np.nan, False + + if multiplicity == 1: + stereo_fov_lon = fov_lon_values[0][1] + stereo_fov_lat = fov_lat_values[0][1] + return stereo_fov_lon, stereo_fov_lat, True + + if ( + multiplicity == 2 + and self.min_ang_diff is not None + and not check_ang_diff(self.min_ang_diff, hillas_psis[0], hillas_psis[1]) + ): + return np.nan, np.nan, False + + n_tel_combs = min(self.n_tel_combinations, multiplicity) + + best_tels_idx = np.arange(multiplicity) + if self.n_best_tels is not None and multiplicity > self.n_best_tels: + best_tels_idx = arg_n_largest(self.n_best_tels, np.array(weights)) + + index_tel_combs = get_combinations(len(best_tels_idx), n_tel_combs) + fov_lons, fov_lats, comb_weights = calc_combs_min_distances( + index_tel_combs, + np.array(fov_lon_values)[best_tels_idx], + np.array(fov_lat_values)[best_tels_idx], + np.array(weights)[best_tels_idx], + ) + + stereo_fov_lon = np.average(fov_lons, weights=comb_weights) + stereo_fov_lat = np.average(fov_lats, weights=comb_weights) + return stereo_fov_lon, stereo_fov_lat, True + + def predict_table(self, mono_predictions: Table) -> Table: + """ + Compute stereo DISP-based direction reconstruction for a full table + of mono predictions. + + This is the table-wise / batch version of the DISP stereo combination. + Each row in ``mono_predictions`` corresponds to a telescope-event. + The function groups rows by subarray event (obs_id, event_id), performs + DISP-based direction reconstruction for each subarray event, and returns + one row per subarray event. + + See :meth:`StereoCombiner.predict_table` for the general + input/output conventions. + """ + prefix_tel = f"{self.prefix}_tel" + + valid = mono_predictions[f"{prefix_tel}_is_valid"].copy() + self._require_disp_column(mono_predictions, prefix_tel) + + obs_ids, event_ids, _, tel_to_array_indices = get_subarray_index( + mono_predictions + ) + n_array_events = len(obs_ids) + + valid = self._apply_min_ang_diff_cut( + mono_predictions=mono_predictions, + valid=valid, + tel_to_array_indices=tel_to_array_indices, + ) + valid = self._apply_n_best_tels_cut( + mono_predictions=mono_predictions, + valid=valid, + tel_to_array_indices=tel_to_array_indices, + ) + + stereo_table = self._init_stereo_table(mono_predictions, obs_ids, event_ids) + alt, az = self._compute_altaz_for_valid( + mono_predictions=mono_predictions, + valid=valid, + tel_to_array_indices=tel_to_array_indices, + n_array_events=n_array_events, + prefix_tel=prefix_tel, + ) + + stereo_table[f"{self.prefix}_alt"] = alt + stereo_table[f"{self.prefix}_az"] = az + stereo_table[f"{self.prefix}_is_valid"] = np.logical_and( + np.isfinite(stereo_table[f"{self.prefix}_alt"]), + np.isfinite(stereo_table[f"{self.prefix}_az"]), + ) + stereo_table[f"{self.prefix}_goodness_of_fit"] = np.full( + len(stereo_table), np.nan + ) + + stereo_table[f"{self.prefix}_telescopes"] = self._collect_telescopes_per_event( + mono_predictions=mono_predictions, + valid=valid, + tel_to_array_indices=tel_to_array_indices, + n_array_events=n_array_events, + ) + + add_defaults_and_meta(stereo_table, _containers[self.property], self.prefix) + return stereo_table + + def _require_disp_column(self, mono_predictions: Table, prefix_tel: str) -> None: + disp_col = f"{prefix_tel}_parameter" + if disp_col not in mono_predictions.colnames: + raise KeyError( + f"Required DISP column '{disp_col}' not found in mono prediction table. " + f"Make sure the mono events were reconstructed with the corresponding " + f"DispReconstructor (prefix='{self.prefix}') before running " + f"StereoDispCombiner." + ) + + def _apply_min_ang_diff_cut( + self, + mono_predictions: Table, + valid: np.ndarray, + tel_to_array_indices: np.ndarray, + ) -> np.ndarray: + if self.min_ang_diff is None: + return valid + + # Reject events with multiplicity 2 and nearly parallel main shower axes + mask_multi2_tels = valid_tels_of_multi(2, tel_to_array_indices[valid]) + if not np.any(mask_multi2_tels): + return valid + + valid_idx = np.flatnonzero(valid) + pairs_in_valid = np.flatnonzero(mask_multi2_tels).reshape(-1, 2) + valid_psis = mono_predictions["hillas_psi"][valid] + + keep_pairs = check_ang_diff( + self.min_ang_diff, + valid_psis[pairs_in_valid[:, 0]], + valid_psis[pairs_in_valid[:, 1]], + ) + if np.all(keep_pairs): + return valid + + tels_to_invalidate = valid_idx[pairs_in_valid[~keep_pairs].ravel()] + valid[tels_to_invalidate] = False + return valid + + def _apply_n_best_tels_cut( + self, + *, + mono_predictions: Table, + valid: np.ndarray, + tel_to_array_indices: np.ndarray, + ) -> np.ndarray: + if self.n_best_tels is None: + return valid + + weights = self._calculate_weights(mono_predictions[valid]) + valid_idx = np.flatnonzero(valid) + valid_tel_to_array_indices = tel_to_array_indices[valid] + + if valid_tel_to_array_indices.size == 0: + return valid + + order = np.lexsort((-np.asarray(weights), valid_tel_to_array_indices)) + + starts = np.empty(valid_tel_to_array_indices.size, dtype=bool) + starts[0] = True + starts[1:] = valid_tel_to_array_indices[1:] != valid_tel_to_array_indices[:-1] + + start_pos = np.where(starts, np.arange(valid_tel_to_array_indices.size), 0) + group_start = np.maximum.accumulate(start_pos) + rank = np.arange(valid_tel_to_array_indices.size) - group_start + + drop_n = rank >= self.n_best_tels + if not np.any(drop_n): + return valid + + valid[valid_idx[order[drop_n]]] = False + return valid + + def _init_stereo_table( + self, mono_predictions: Table, obs_ids: np.ndarray, event_ids: np.ndarray + ) -> Table: + stereo_table = Table({"obs_id": obs_ids, "event_id": event_ids}) + for colname in ("obs_id", "event_id"): + stereo_table[colname].description = mono_predictions[colname].description + return stereo_table + + def _compute_altaz_for_valid( + self, + mono_predictions: Table, + valid: np.ndarray, + tel_to_array_indices: np.ndarray, + n_array_events: int, + prefix_tel: str, + ): + if np.count_nonzero(valid) == 0: + nan = u.Quantity( + np.full(n_array_events, np.nan), u.deg, copy=COPY_IF_NEEDED + ) + return nan, nan + + weights = self._calculate_weights(mono_predictions[valid]) + _, _, valid_multiplicity, _ = get_subarray_index(mono_predictions[valid]) + + fov_lon_values, fov_lat_values = calc_fov_lon_lat( + mono_predictions[valid], prefix_tel + ) + + combs_array, combs_to_multi_indices = create_combs_array( + valid_multiplicity.max(), self.n_tel_combinations + ) + index_tel_combs, n_combs = get_index_combs( + valid_multiplicity, + combs_array, + combs_to_multi_indices, + self.n_tel_combinations, + ) + combs_to_array_indices = np.repeat(np.arange(len(valid_multiplicity)), n_combs) + + comb_fov_lons, comb_fov_lats, comb_weights = calc_combs_min_distances( + index_tel_combs, + fov_lon_values, + fov_lat_values, + weights, + ) + + all_valid = np.ones(len(comb_weights), dtype=bool) + fov_lon_combs_mean, _ = weighted_mean_std_ufunc( + comb_fov_lons, + all_valid, + combs_to_array_indices, + n_combs, + weights=comb_weights, + ) + fov_lat_combs_mean, _ = weighted_mean_std_ufunc( + comb_fov_lats, + all_valid, + combs_to_array_indices, + n_combs, + weights=comb_weights, + ) + + valid_tel_to_array_indices = tel_to_array_indices[valid] + valid_array_indices = np.unique(valid_tel_to_array_indices) + + if self.n_tel_combinations >= 3: + fill_lower_multiplicities( + fov_lon_combs_mean, + fov_lat_combs_mean, + self.n_tel_combinations, + valid_tel_to_array_indices, + valid_multiplicity, + fov_lon_values, + fov_lat_values, + weights, + ) + + fov_lon_mean = np.full(n_array_events, np.nan) + fov_lat_mean = np.full(n_array_events, np.nan) + fov_lon_mean[valid_array_indices] = fov_lon_combs_mean + fov_lat_mean[valid_array_indices] = fov_lat_combs_mean + + _, indices_first_tel_in_array = np.unique( + tel_to_array_indices, return_index=True + ) + alt, az = telescope_to_horizontal( + lon=u.Quantity(fov_lon_mean, u.deg, copy=COPY_IF_NEEDED), + lat=u.Quantity(fov_lat_mean, u.deg, copy=COPY_IF_NEEDED), + pointing_alt=u.Quantity( + mono_predictions["subarray_pointing_lat"][indices_first_tel_in_array], + u.deg, + copy=COPY_IF_NEEDED, + ), + pointing_az=u.Quantity( + mono_predictions["subarray_pointing_lon"][indices_first_tel_in_array], + u.deg, + copy=COPY_IF_NEEDED, + ), + ) + + # Fill single telescope events from mono_predictions + index_single_tel_events = valid_array_indices[valid_multiplicity == 1] + if index_single_tel_events.size > 0: + mask_single_tel_events = valid_tels_of_multi(1, valid_tel_to_array_indices) + alt[index_single_tel_events] = mono_predictions[f"{prefix_tel}_alt"][valid][ + mask_single_tel_events + ] + az[index_single_tel_events] = mono_predictions[f"{prefix_tel}_az"][valid][ + mask_single_tel_events + ] + + return alt, az + + def _collect_telescopes_per_event( + self, + mono_predictions: Table, + valid: np.ndarray, + tel_to_array_indices: np.ndarray, + n_array_events: int, + ): + tel_ids = [[] for _ in range(n_array_events)] + for index, tel_id in zip( + tel_to_array_indices[valid], mono_predictions["tel_id"][valid] + ): + tel_ids[index].append(tel_id) + return tel_ids diff --git a/src/ctapipe/reco/telescope_event_handling.py b/src/ctapipe/reco/telescope_event_handling.py index 5e9da9621e8..9d4dc8db545 100644 --- a/src/ctapipe/reco/telescope_event_handling.py +++ b/src/ctapipe/reco/telescope_event_handling.py @@ -1,12 +1,30 @@ """Helper functions for array-event-wise aggregation of telescope events.""" +from functools import lru_cache +from itertools import combinations, product + +import astropy.units as u import numpy as np from numba import njit, uint64 -__all__ = ["get_subarray_index", "weighted_mean_std_ufunc"] +from ..core.env import CTAPIPE_DISABLE_NUMBA_CACHE + +__all__ = [ + "get_subarray_index", + "weighted_mean_std_ufunc", + "get_combinations", + "binary_combinations", + "check_ang_diff", + "calc_combs_min_distances", + "valid_tels_of_multi", + "fill_lower_multiplicities", + "calc_fov_lon_lat", + "create_combs_array", + "get_index_combs", +] -@njit +@njit(cache=not CTAPIPE_DISABLE_NUMBA_CACHE) def _get_subarray_index(obs_ids, event_ids): n_tel_events = len(obs_ids) idx = np.zeros(n_tel_events, dtype=uint64) @@ -48,26 +66,34 @@ def _get_subarray_index(obs_ids, event_ids): def get_subarray_index(tel_table): """ - Get the subarray-event-wise information from a table of telescope events. + Get 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. + Extract the unique subarray events contained in a table of telescope-event + rows and return their observation IDs, event IDs, multiplicities, and an + index mapping from each telescope-event row to its corresponding subarray + event. - This requires the telescope events of one subarray event to be come - in a single block. + This requires that all telescope events belonging to the same subarray + event appear in a single contiguous block in ``tel_table``. Parameters ---------- - tel_table: astropy.table.Table - table with telescope events as rows + tel_table : astropy.table.Table + Table with telescope events as rows. Must contain the columns + ``obs_id`` and ``event_id`` and be ordered such that telescope rows + of the same (obs_id, event_id) are contiguous. 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_ids : np.ndarray + Observation IDs of the subarray events. One entry per subarray event. + event_ids : np.ndarray + Event IDs of the subarray events. One entry per subarray event. + multiplicity : np.ndarray + Number of telescope-event rows contributing to each subarray event. + tel_to_array_index : np.ndarray + Integer array mapping each telescope-event row to the corresponding + subarray-event index. Has the same length as ``tel_table``. """ obs_idx = tel_table["obs_id"] event_idx = tel_table["event_id"] @@ -76,9 +102,26 @@ def get_subarray_index(tel_table): 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. + Sum telescope-event values per subarray event. + + Groups telescope-event values by their corresponding subarray-event index + and computes the group-wise sum using ``np.add.at``. + + Parameters + ---------- + tel_data : np.ndarray + Values for each telescope event (one value per telescope-event row). + n_array_events : int + Total number of subarray events (size of the grouped output). + indices : np.ndarray + Integer array mapping each telescope-event row to its corresponding + subarray-event index. Must have the same length as ``tel_data``. + + Returns + ------- + summed : np.ndarray + Array of shape ``(n_array_events,)`` containing the sum of ``tel_data`` + over all telescope-event rows belonging to each subarray event. """ combined_values = np.zeros(n_array_events) np.add.at(combined_values, indices, tel_data) @@ -93,28 +136,44 @@ def weighted_mean_std_ufunc( weights=np.array([1]), ): """ - Calculate the weighted mean and standard deviation for each array event - over the corresponding telescope events. + Compute weighted mean and standard deviation per subarray event. + + Telescope-event values (``tel_values``) are grouped by subarray event using + ``indices`` (mapping each telescope-event row to an subarray-event index). + For each subarray event, the function computes the weighted mean and the + weighted standard deviation over the corresponding telescope-event rows. + + Invalid telescope rows are excluded using ``valid_tel``. Subarray events + without any valid telescope contribution return ``NaN`` for both mean and + standard deviation. Parameters ---------- - tel_values: np.ndarray - values for each telescope event - valid_tel: array-like - boolean mask giving the valid values of ``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) + tel_values : np.ndarray + Values for each telescope event (one value per telescope-event row). + valid_tel : array-like + Boolean mask selecting valid telescope-event rows in ``tel_values``. + Must have the same length as ``tel_values``. + indices : np.ndarray + Integer array mapping each telescope-event row to the corresponding + subarray-event index. This is the fourth return value of + ``get_subarray_index``. Must have the same length as ``tel_values``. + multiplicity : np.ndarray + Number of telescope-event rows per subarray event (one entry per + subarray event), in the same order as the subarray events encoded in + ``indices``. + weights : np.ndarray, optional + Weights used for averaging. Must be broadcastable to ``tel_values``. + If not provided, equal weights are used. Returns ------- - Tuple(np.ndarray, np.ndarray) - weighted mean and standard deviation for each array event + mean : np.ndarray + Weighted mean value for each subarray event. Entries are + ``NaN`` for events without valid telescope contributions. + std : np.ndarray + Weighted standard deviation for each subarray event. + Entries are ``NaN`` where the mean is undefined. """ n_array_events = len(multiplicity) # avoid numerical problems by very large or small weights @@ -136,11 +195,490 @@ def weighted_mean_std_ufunc( valid = sum_of_weights > 0 mean[valid] = sum_prediction[valid] / sum_of_weights[valid] - sum_sq_residulas = _grouped_add( + sum_sq_residuals = _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] + variance[valid] = sum_sq_residuals[valid] / sum_of_weights[valid] return mean, np.sqrt(variance) + + +@lru_cache(maxsize=4096) +def get_combinations(array_length: int, comb_size: int) -> np.ndarray: + """ + Generate all index combinations of a fixed size. + + Returns all ``comb_size``-element combinations chosen from the index range + ``0 .. array_length - 1``. Results are cached to speed up repeated calls + with the same arguments. + + Parameters + ---------- + array_length : int + Length of the index range to draw from. + comb_size : int + Size of each combination. + + Returns + ------- + combs : np.ndarray + Integer array of shape ``(n_combs, comb_size)`` containing all index + combinations, where ``n_combs = binom(array_length, comb_size)``. + """ + return np.array(list(combinations(range(array_length), comb_size))) + + +@lru_cache(maxsize=4096) +def binary_combinations(k: int) -> np.ndarray: + """ + Generate all binary (0/1) vectors of length ``k``. + + Returns the full set of binary sign assignments used e.g. to enumerate + DISP head–tail sign choices. Results are cached to speed up repeated calls. + + Parameters + ---------- + k : int + Length of each binary vector. + + Returns + ------- + combs : np.ndarray + Integer array of shape ``(2**k, k)`` containing all possible binary + combinations. Each row is one assignment. + """ + return np.array(list(product([0, 1], repeat=k)), dtype=int) + + +def check_ang_diff(min_ang_diff, psi1, psi2): + """ + Check whether two Hillas main-axis orientations are sufficiently separated. + + Computes an axis-invariant angular difference between two orientation angles, + treating directions that differ by 180° as parallel. The resulting difference + is folded into the interval [0°, 90°] and compared against ``min_ang_diff``. + + Parameters + ---------- + min_ang_diff : float + Minimum required angular separation in degrees. + psi1, psi2 : astropy.units.Quantity + Orientation angles of the Hillas main axes. Scalars or array-like objects + of identical shape with angular units (e.g. degrees). + + Returns + ------- + keep : bool or np.ndarray + Boolean value or boolean array indicating whether the (axis-invariant) + angular separation is greater than or equal to ``min_ang_diff``. + ``True`` means the event should be kept; ``False`` means the axes are + too parallel. + """ + ang_diff = np.abs(psi1 - psi2) % (180 * u.deg) + ang_diff = np.minimum(ang_diff, 180 * u.deg - ang_diff) + return ang_diff >= (min_ang_diff * u.deg) + + +def calc_combs_min_distances(index_tel_combs, fov_lon_values, fov_lat_values, weights): + """ + Select the DISP sign assignment that minimizes the weighted SSE per combination. + + For each telescope combination (size ``k``), evaluate all ``2**k`` possible + binary sign assignments (DISP head–tail ambiguity). For each assignment, + compute the weighted mean FoV direction and the weighted sum of squared + distances (SSE) of the telescope directions to this mean. The assignment + with minimal SSE is chosen. + + The function returns the weighted mean FoV longitude/latitude for the best + assignment of each telescope combination, plus the summed weight per + combination. + + Parameters + ---------- + index_tel_combs : np.ndarray + Integer array of shape ``(n_combs, k)`` containing index combinations of + telescope events forming each combination. + fov_lon_values : np.ndarray + Array of shape ``(n_tel_events, 2)`` containing the two possible FoV + longitude values (SIGN = 0/1, i.e. DISP sign choices) for each telescope + event, in degrees. + fov_lat_values : np.ndarray + Array of shape ``(n_tel_events, 2)`` containing the two possible FoV + latitude values (SIGN = 0/1) for each telescope event, in degrees. + weights : np.ndarray + Array of shape ``(n_tel_events,)`` containing the weight per telescope + event. + + Returns + ------- + weighted_lons : np.ndarray + Array of shape ``(n_combs,)`` with the weighted mean FoV longitudes of + the best sign assignment for each telescope combination. + weighted_lats : np.ndarray + Array of shape ``(n_combs,)`` with the weighted mean FoV latitudes of + the best sign assignment for each telescope combination. + combined_weights : np.ndarray + Array of shape ``(n_combs,)`` containing the sum of telescope weights + contributing to each combination. + """ + mapped_weights = weights[index_tel_combs] + combined_weights = mapped_weights.sum(axis=1) + + _, k = index_tel_combs.shape + sign_combs = binary_combinations(k) + + lon_combs = fov_lon_values[index_tel_combs] # (n_combs, k, 2) + lat_combs = fov_lat_values[index_tel_combs] # (n_combs, k, 2) + + lon0 = lon_combs[..., 0] + lon1 = lon_combs[..., 1] + lat0 = lat_combs[..., 0] + lat1 = lat_combs[..., 1] + + s = sign_combs.T[None, :, :] # (1, k, 2**k), 0/1 + lon_vals = lon0[..., None] * (1 - s) + lon1[..., None] * s + lat_vals = lat0[..., None] * (1 - s) + lat1[..., None] * s + + w = mapped_weights[:, :, None] + wsum = w.sum(axis=1) + + lon_mu = (lon_vals * w).sum(axis=1) / wsum + lat_mu = (lat_vals * w).sum(axis=1) / wsum + + dlon = lon_vals - lon_mu[:, None, :] + dlat = lat_vals - lat_mu[:, None, :] + sse = (w * (dlon * dlon + dlat * dlat)).sum(axis=1) + + argmin = np.argmin(sse, axis=1) + best_signs = sign_combs[argmin] + + best_lon = lon0 * (1 - best_signs) + lon1 * best_signs + best_lat = lat0 * (1 - best_signs) + lat1 * best_signs + + weighted_lons = np.average(best_lon, weights=mapped_weights, axis=1) + weighted_lats = np.average(best_lat, weights=mapped_weights, axis=1) + + return weighted_lons, weighted_lats, combined_weights + + +def valid_tels_of_multi(multi, valid_tel_to_array_indices): + """ + Create a boolean mask selecting telescope-event rows that belong to + subarray events with a given multiplicity. + + The function assumes that telescope events are grouped contiguously + by subarray event, i.e. all rows belonging to the same subarray event appear + in a single consecutive block in ``valid_tel_to_array_indices``. + Under this assumption, the multiplicity of an subarray event corresponds + to the length of its contiguous block. + + For all subarray events whose block length equals ``multi``, the function + returns a boolean mask that is ``True`` for all telescope rows belonging + to those events and ``False`` otherwise. + + This block-based approach avoids repeated membership tests (e.g. + ``np.isin``) and is faster for large arrays. + + Parameters + ---------- + multi : int + Target multiplicity (number of telescope events per subarray event). + valid_tel_to_array_indices : np.ndarray + One-dimensional array mapping each valid telescope event to its + corresponding subarray-event index. Must be ordered such that + telescope events of the same subarray event are contiguous. + + Returns + ------- + mask : np.ndarray + Boolean array of the same length as ``valid_tel_to_array_indices``. + Entries are ``True`` for telescope-event rows belonging to subarray + events with multiplicity ``multi`` and ``False`` otherwise. + + """ + change = np.empty(len(valid_tel_to_array_indices), dtype=bool) + change[0] = True + change[1:] = valid_tel_to_array_indices[1:] != valid_tel_to_array_indices[:-1] + + starts = np.flatnonzero(change) + lengths = np.diff(np.append(starts, len(valid_tel_to_array_indices))) + good_blocks = lengths == multi + mask = np.repeat(good_blocks, lengths) + + return mask + + +def fill_lower_multiplicities( + fov_lon_combs_mean, + fov_lat_combs_mean, + n_tel_combinations, + valid_tel_to_array_indices, + valid_multiplicity, + fov_lon_values, + fov_lat_values, + weights, +): + """ + Fill stereo FoV estimates for events with multiplicity < ``n_tel_combinations``. + + For subarray events whose multiplicity is smaller than ``n_tel_combinations`` + but at least 2, recompute the stereo direction using *all* available + telescopes of the event (i.e. combination size equals the event multiplicity). + The optimal DISP sign assignment is determined via + :func:`calc_combs_min_distances`. + + Results are written in-place into ``fov_lon_combs_mean`` and + ``fov_lat_combs_mean`` at the positions of the affected events. + + Parameters + ---------- + fov_lon_combs_mean : np.ndarray + Array of shape ``(n_array_events,)`` holding the mean FoV longitudes per + subarray event. Updated in-place for lower-multiplicity events. + fov_lat_combs_mean : np.ndarray + Array of shape ``(n_array_events,)`` holding the mean FoV latitudes per + subarray event. Updated in-place for lower-multiplicity events. + n_tel_combinations : int + Nominal number of telescopes used per stereo combination. + valid_tel_to_array_indices : np.ndarray + One-dimensional array mapping each valid telescope-event row to its + corresponding subarray-event index. Telescope rows must be grouped + contiguously by subarray event. + valid_multiplicity : np.ndarray + Multiplicity per valid subarray event (one entry per subarray event), aligned + with the indices used in ``valid_tel_to_array_indices``. + fov_lon_values : np.ndarray + Array of shape ``(n_valid_tel_events, 2)`` containing the two possible + FoV longitude values (DISP sign choices) for each valid telescope event. + fov_lat_values : np.ndarray + Array of shape ``(n_valid_tel_events, 2)`` containing the two possible + FoV latitude values (DISP sign choices) for each valid telescope event. + weights : np.ndarray + Array of shape ``(n_valid_tel_events,)`` with the weight per valid + telescope event. + + Returns + ------- + None + Operates in-place and does not return a value. + """ + for multi in range(n_tel_combinations - 1, 1, -1): + multi_mask = valid_multiplicity == multi + if not np.any(multi_mask): + continue + + mask = valid_tels_of_multi(multi, valid_tel_to_array_indices) + index_tel_combs = np.arange(len(valid_tel_to_array_indices))[mask].reshape( + -1, multi + ) + + lons, lats, _ = calc_combs_min_distances( + index_tel_combs, + fov_lon_values, + fov_lat_values, + weights, + ) + fov_lon_combs_mean[multi_mask] = lons + fov_lat_combs_mean[multi_mask] = lats + + +def calc_fov_lon_lat(tel_table, prefix="DispReconstructor_tel"): + """ + Compute the two possible FoV longitude/latitude coordinates per telescope event. + + For each telescope event, compute two candidate shower directions in the + telescope FoV corresponding to the two possible DISP sign assignments + (SIGN = −1 and +1). The candidates are obtained by shifting the Hillas + centroid along the main shower axis by the absolute DISP distance. + + Parameters + ---------- + tel_table : astropy.table.Table + Table containing Hillas parameters and DISP reconstruction results. + Must include: + - ``hillas_fov_lon`` : Hillas centroid longitude (Quantity) + - ``hillas_fov_lat`` : Hillas centroid latitude (Quantity) + - ``hillas_psi`` : Hillas orientation angle (Quantity) + - ``_parameter`` : DISP distance from the centroid (float) + prefix : str, optional + Prefix used to access the DISP parameter column + (default: "DispReconstructor_tel"). + + Returns + ------- + fov_lon : np.ndarray + Array of shape ``(n_tel_events, 2)`` containing the candidate FoV + longitudes (in degrees) for each telescope event (DISP signs −1 and +1). + fov_lat : np.ndarray + Array of shape ``(n_tel_events, 2)`` containing the candidate FoV + latitudes (in degrees) for each telescope event (DISP signs −1 and +1). + """ + hillas_fov_lon = tel_table["hillas_fov_lon"].quantity.to_value(u.deg) + hillas_fov_lat = tel_table["hillas_fov_lat"].quantity.to_value(u.deg) + hillas_psi = tel_table["hillas_psi"].quantity.to_value(u.rad) + disp = tel_table[f"{prefix}_parameter"] + signs = np.array([-1, 1]) + + cos_psi = np.cos(hillas_psi) + sin_psi = np.sin(hillas_psi) + abs_disp = np.abs(disp)[:, None] + lons = hillas_fov_lon[:, None] + signs * abs_disp * cos_psi[:, None] + lats = hillas_fov_lat[:, None] + signs * abs_disp * sin_psi[:, None] + + return lons, lats + + +def create_combs_array(max_multiplicity, k): + """ + Precompute all possible ``k``-combinations for multiplicities up to + ``max_multiplicity``. + + For each multiplicity ``m`` in the range ``k .. max_multiplicity``, this + function generates all index combinations of size ``k`` chosen from + ``0 .. m-1`` (using ``get_combinations(m, k)``) and concatenates them into + a single array. + + In addition, it returns an array mapping each combination row to the + multiplicity ``m`` it was generated from. This mapping can later be used + to quickly select the correct subset of combinations for a given event + multiplicity without recomputing combinations. + + Parameters + ---------- + max_multiplicity : int + Maximum multiplicity to consider (inclusive). + k : int + Size of the combinations (must satisfy ``k >= 2`` and + ``k <= max_multiplicity``). + + Returns + ------- + combs_array : np.ndarray + Array of shape ``(n_total_combinations, k)`` containing all possible + ``k``-combinations for multiplicities ranging from ``k`` up to + ``max_multiplicity``. Combinations are ordered by increasing + multiplicity. + combs_to_multi_indices : np.ndarray + One-dimensional integer array mapping each row in ``combs_array`` to + the multiplicity it was generated from. Has the same length as + ``combs_array``. + """ + combs_array = get_combinations(k, k) + for i in range(k + 1, max_multiplicity + 1): + combs = get_combinations(i, k) + combs_array = np.concatenate([combs_array, combs]) + + n_combs = _calc_n_combs(np.arange(k, max_multiplicity + 1), k) + combs_to_multi_indices = np.repeat(np.arange(k, max_multiplicity + 1), n_combs) + + return combs_array, combs_to_multi_indices + + +@njit(cache=not CTAPIPE_DISABLE_NUMBA_CACHE) +def _binomial(n, k): + """ + Compute the binomial coefficient (`n` choose `k`). + + Parameters + ---------- + n : int + Total number of items. + k : int + Number of selections. + + Returns + ------- + int + The binomial coefficient (n choose k). + """ + if k > n or k < 0: + return 0 + k = min(k, n - k) + c = 1 + for i in range(k): + c = c * (n - i) // (i + 1) + return c + + +@njit(cache=not CTAPIPE_DISABLE_NUMBA_CACHE) +def _calc_n_combs(multiplicity, k): + """ + Compute the number of ``k``-combinations for each multiplicity. + + Parameters + ---------- + multiplicity : np.ndarray + One-dimensional array of multiplicities (number of telescope-event rows) + for each subarray event. + k : int + Size of the combinations. + + Returns + ------- + n_combs : np.ndarray + One-dimensional integer array containing the number of possible + ``k``-combinations for each entry in ``multiplicity``. + """ + n_combs = np.empty(len(multiplicity), dtype=np.int64) + for i in range(len(multiplicity)): + n_combs[i] = _binomial(multiplicity[i], k) + + return n_combs + + +@njit(cache=not CTAPIPE_DISABLE_NUMBA_CACHE) +def get_index_combs(multiplicities, combs_array, combs_to_multi_indices, k): + """ + Build telescope-event index combinations for all subarray events. + + For each subarray event with multiplicity ``m = multiplicities[i]``, this + function selects the precomputed ``k``-combinations corresponding to ``m`` + from ``combs_array`` (using ``combs_to_multi_indices``) and offsets them + into the global telescope-event indexing scheme. + + The global telescope-event ordering is assumed to be the concatenation of + telescope-event blocks per subarray event. Under this assumption, the start + index of event ``i`` equals the cumulative sum of previous multiplicities, + and combinations for event ``i`` can be obtained by adding that offset. + + Parameters + ---------- + multiplicities : np.ndarray + One-dimensional array of multiplicities (number of telescope-event rows) + for each subarray event. + combs_array : np.ndarray + Precomputed combination indices for multiplicities in the range + ``k .. max_multiplicity``. Usually produced by ``create_combs_array``. + Shape: ``(n_total_combinations, k)``. + combs_to_multi_indices : np.ndarray + One-dimensional array mapping each row in ``combs_array`` to the + multiplicity it was generated from. Same length as ``combs_array``. + k : int + Size of the combinations. + + Returns + ------- + index_tel_combs : np.ndarray + Array of shape ``(n_total_combinations_over_events, k)`` containing the + telescope-event indices for all ``k``-combinations across all subarray + events. Indices refer to the flattened telescope-event ordering. + n_combs : np.ndarray + One-dimensional array giving the number of ``k``-combinations for each + subarray event, in the same order as ``multiplicities``. + """ + n_combs = _calc_n_combs(multiplicities, k) + total_combs = np.sum(n_combs) + index_tel_combs = np.empty((total_combs, k), dtype=np.int64) + cum_multiplicities = 0 + idx = 0 + + for i in range(len(multiplicities)): + mask = combs_to_multi_indices == multiplicities[i] + selected_combs = combs_array[mask] + cum_multiplicities + index_tel_combs[idx : idx + len(selected_combs)] = selected_combs + idx += len(selected_combs) + cum_multiplicities += multiplicities[i] + + return index_tel_combs, n_combs diff --git a/src/ctapipe/reco/tests/test_stereo_combination.py b/src/ctapipe/reco/tests/test_stereo_combination.py index 59127337848..8fff9868425 100644 --- a/src/ctapipe/reco/tests/test_stereo_combination.py +++ b/src/ctapipe/reco/tests/test_stereo_combination.py @@ -6,15 +6,20 @@ from ctapipe.containers import ( ArrayEventContainer, + ArrayPointingContainer, + DispContainer, HillasParametersContainer, ImageParametersContainer, + MorphologyContainer, ParticleClassificationContainer, ReconstructedContainer, ReconstructedEnergyContainer, ReconstructedGeometryContainer, + TelescopeReconstructedContainer, ) +from ctapipe.core.traits import TraitError from ctapipe.reco.reconstructor import ReconstructionProperty -from ctapipe.reco.stereo_combination import StereoMeanCombiner +from ctapipe.reco.stereo_combination import StereoDispCombiner, StereoMeanCombiner @pytest.fixture(scope="module") @@ -25,13 +30,18 @@ def mono_table(): """ return Table( { - "obs_id": [1, 1, 1, 1, 1, 2], - "event_id": [1, 1, 1, 2, 2, 1], - "tel_id": [1, 2, 3, 5, 7, 1], - "hillas_intensity": [1, 2, 0, 1, 5, 9], - "hillas_width": [0.1, 0.2, 0.1, 0.1, 0.2, 0.1] * u.deg, - "hillas_length": 3 * ([0.1, 0.2, 0.1, 0.1, 0.2, 0.1] * u.deg), - "dummy_tel_energy": [1, 10, 4, 0.5, 0.7, 1] * u.TeV, + "obs_id": [1, 1, 1, 1, 1, 2, 2, 2, 2, 2], + "event_id": [1, 1, 1, 2, 2, 1, 2, 2, 2, 2], + "tel_id": [1, 2, 3, 5, 7, 1, 1, 3, 4, 5], + "hillas_intensity": [1, 2, 0, 1, 5, 9, 10, 20, 1, 2], + "hillas_width": [0.1, 0.2, 0.1, 0.1, 0.2, 0.1, 0.1, 0.2, 0.1, 0.2] * u.deg, + "hillas_length": 3 + * ([0.1, 0.2, 0.1, 0.1, 0.2, 0.1, 0.1, 0.2, 0.1, 0.2] * u.deg), + "hillas_fov_lon": [-0.5, 0, 0.5, -1, 1, 1.5, -0.5, 0, 0.5, -1] * u.deg, + "hillas_fov_lat": [0.3, -0.3, 0.3, 0.5, 0.5, 0.2, 0.3, -0.3, 0.3, 0.5] + * u.deg, + "hillas_psi": [40, 85, -40, -35, 35, 55, 40, 85, -40, -35] * u.deg, + "dummy_tel_energy": [1, 10, 4, 0.5, 0.7, 1, 1, 9, 4, 0.5] * u.TeV, "dummy_tel_is_valid": [ True, True, @@ -39,8 +49,12 @@ def mono_table(): True, False, False, + True, + True, + True, + True, ], - "classifier_tel_prediction": [1, 0, 0.5, 0, 0.6, 1], + "classifier_tel_prediction": [1, 0, 0.5, 0, 0.6, 1, 1, 0, 0.5, 0], "classifier_tel_is_valid": [ True, True, @@ -48,17 +62,41 @@ def mono_table(): True, True, True, + True, + True, + True, + True, ], - "disp_tel_alt": [58.5, 58, 62.5, 72, 74.5, 81] * u.deg, - "disp_tel_az": [12.5, 15, 13, 21, 20, 14.5] * u.deg, + "disp_tel_alt": [58.5, 58, 62.5, 72, 74.5, 81, 58.5, 58, 62.5, 72] * u.deg, + "disp_tel_az": [12.5, 15, 13, 21, 20, 14.5, 12.5, 15, 13, 21] * u.deg, "disp_tel_is_valid": [ True, False, True, + False, + False, True, True, True, + True, + True, + ], + "disp_tel_parameter": [0.65, 1.1, 0.7, 0.9, 1, 0.5, 0.65, 1.1, 0.7, 0.9] + * u.deg, + "disp_tel_sign_score": [ + 0.65, + 0.87, + 0.7, + 0.5, + 0.3, + 0.99, + 0.65, + 0.3, + 0.2, + 0.86, ], + "subarray_pointing_lat": 10 * [70] * u.deg, + "subarray_pointing_lon": 10 * [0] * u.deg, } ) @@ -80,25 +118,32 @@ def test_predict_mean_energy(weights, mono_table): "dummy_goodness_of_fit", "dummy_telescopes", ] - assert_array_equal(stereo["obs_id"], np.array([1, 1, 2])) - assert_array_equal(stereo["event_id"], np.array([1, 2, 1])) + assert_array_equal(stereo["obs_id"], np.array([1, 1, 2, 2])) + assert_array_equal(stereo["event_id"], np.array([1, 2, 1, 2])) if weights == "intensity": - assert_array_equal(stereo["dummy_energy"], [7, 0.5, np.nan] * u.TeV) + assert_allclose( + stereo["dummy_energy"].quantity, + [7, 0.5, np.nan, 5.90909091] * u.TeV, + atol=1e-7, + ) assert_allclose( stereo["dummy_energy_uncert"].quantity, - [4.242641, 0, np.nan] * u.TeV, + [4.242641, 0, np.nan, 3.869959] * 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"].quantity, [5, 0.5, np.nan, 3.625] * u.TeV, atol=1e-7 + ) assert_allclose( stereo["dummy_energy_uncert"].quantity, - [3.741657, 0, np.nan] * u.TeV, + [3.741657, 0, np.nan, 3.3796265] * 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) + assert_array_equal(stereo["dummy_telescopes"][3], np.array([1, 3, 4, 5])) def test_predict_mean_classification(mono_table): @@ -115,16 +160,17 @@ def test_predict_mean_classification(mono_table): "classifier_goodness_of_fit", "classifier_telescopes", ] - assert_array_equal(stereo["obs_id"], np.array([1, 1, 2])) - assert_array_equal(stereo["event_id"], np.array([1, 2, 1])) + assert_array_equal(stereo["obs_id"], np.array([1, 1, 2, 2])) + assert_array_equal(stereo["event_id"], np.array([1, 2, 1, 2])) assert_array_equal( stereo["classifier_prediction"], - [0.5, 0.3, 1], + [0.5, 0.3, 1, 0.375], ) tel_ids = stereo["classifier_telescopes"] assert_array_equal(tel_ids[0], [1, 2]) assert_array_equal(tel_ids[1], [5, 7]) assert_array_equal(tel_ids[2], [1]) + assert_array_equal(tel_ids[3], [1, 3, 4, 5]) def test_predict_mean_disp(mono_table): @@ -142,22 +188,23 @@ def test_predict_mean_disp(mono_table): assert "obs_id" in stereo.colnames assert "event_id" in stereo.colnames - assert_array_equal(stereo["obs_id"], np.array([1, 1, 2])) - assert_array_equal(stereo["event_id"], np.array([1, 2, 1])) + assert_array_equal(stereo["obs_id"], np.array([1, 1, 2, 2])) + assert_array_equal(stereo["event_id"], np.array([1, 2, 1, 2])) assert_allclose( stereo["disp_alt"].quantity, - [60.5002328, 73.2505989, 81] * u.deg, + [60.5002328, np.nan, 81, 62.773741] * u.deg, atol=1e-7, ) assert_allclose( stereo["disp_az"].quantity, - [12.7345693, 20.5362510, 14.5] * u.deg, + [12.7345693, np.nan, 14.5, 14.792156] * u.deg, atol=1e-7, ) tel_ids = stereo["disp_telescopes"] assert_array_equal(tel_ids[0], [1, 3]) - assert_array_equal(tel_ids[1], [5, 7]) + assert_array_equal(tel_ids[1], []) assert_array_equal(tel_ids[2], [1]) + assert_array_equal(tel_ids[3], [1, 3, 4, 5]) @pytest.mark.parametrize("weights", ["aspect-weighted-intensity", "intensity", "none"]) @@ -170,7 +217,10 @@ def test_mean_prediction_single_event(weights): intensity=intensity, width=0.1 * u.deg, length=0.3 * u.deg, - ) + ), + morphology=MorphologyContainer( + n_pixels=10, + ), ) event.dl2.tel[25] = ReconstructedContainer( @@ -252,3 +302,345 @@ def test_reconstructed_container_warning(): with pytest.warns(CTAPipeDeprecationWarning, match="renamed"): container.classification = ParticleClassificationContainer() + + +def _make_disp_event(event_dict, prefix="dummy"): + event = ArrayEventContainer() + + for i in range(len(event_dict["tel_id"])): + event.dl1.tel[event_dict["tel_id"][i]].parameters = ImageParametersContainer( + hillas=HillasParametersContainer( + intensity=event_dict["hillas_intensity"][i], + fov_lon=event_dict["hillas_fov_lon"][i], + fov_lat=event_dict["hillas_fov_lat"][i], + psi=event_dict["hillas_psi"][i], + width=event_dict["hillas_width"][i], + length=event_dict["hillas_length"][i], + ), + morphology=MorphologyContainer( + n_pixels=10, + ), + ) + + event.dl2.tel[event_dict["tel_id"][i]] = TelescopeReconstructedContainer( + disp={ + prefix: DispContainer( + parameter=event_dict["disp_tel_parameter"][i], + sign_score=event_dict["disp_tel_sign_score"][i], + ) + }, + geometry={ + prefix: ReconstructedGeometryContainer( + alt=event_dict["disp_tel_alt"][i], + az=event_dict["disp_tel_az"][i], + is_valid=event_dict["disp_tel_is_valid"][i], + ) + }, + ) + + event.monitoring.pointing = ArrayPointingContainer( + array_azimuth=0 * u.deg, array_altitude=70 * u.deg + ) + return event + + +def test_disp_combiner_trait_validation(): + with pytest.raises(TraitError): + StereoDispCombiner(n_tel_combinations=1) + + with pytest.raises(TraitError): + StereoDispCombiner(n_best_tels=1) + + with pytest.raises(TraitError): + StereoDispCombiner(n_tel_combinations=3, n_best_tels=2) + + +@pytest.mark.parametrize("weights", ["aspect-weighted-intensity", "intensity", "none"]) +def test_disp_combiner_single_event(weights): + event_dict = { + "tel_id": [1, 2, 9, 10], + "hillas_intensity": [100, 200, 75, 30], + "hillas_width": [0.1, 0.2, 0.1, 0.1] * u.deg, + "hillas_length": 3 * ([0.1, 0.2, 0.1, 0.1] * u.deg), + "hillas_fov_lon": [-0.5, 0, 0.5, 0.1] * u.deg, + "hillas_fov_lat": [0.3, -0.3, 0.3, 0.2] * u.deg, + "hillas_psi": [40, 85, -40, 30] * u.deg, + "disp_tel_alt": [58.5, 58, 62.5, 20] * u.deg, + "disp_tel_az": [12.5, 15, 13, 30] * u.deg, + "disp_tel_parameter": [0.65, 1.1, 0.7, 1.0] * u.deg, + "disp_tel_sign_score": [0.95, 0.98, 0.66, 0], + "disp_tel_is_valid": [True, True, True, False], + } + + event = _make_disp_event(event_dict) + + disp_combiner = StereoDispCombiner( + prefix="dummy", + weights=weights, + ) + disp_combiner(event) + if weights in ["intensity", "aspect-weighted-intensity"]: + assert u.isclose(event.dl2.stereo.geometry["dummy"].alt, 70.76579427 * u.deg) + assert u.isclose(event.dl2.stereo.geometry["dummy"].az, 0.13152707 * u.deg) + elif weights == "none": + assert u.isclose(event.dl2.stereo.geometry["dummy"].alt, 70.75451665 * u.deg) + assert u.isclose(event.dl2.stereo.geometry["dummy"].az, 0.05821327 * u.deg) + + +def test_disp_combiner_single_event_disp_parameter(): + event = ArrayEventContainer() + event.dl2.tel[1] = TelescopeReconstructedContainer( + geometry={ + "dummy": ReconstructedGeometryContainer( + is_valid=True, + ) + }, + ) + + disp_combiner = StereoDispCombiner( + prefix="dummy", + ) + + with pytest.raises(RuntimeError, match="No valid DISP reconstruction parameter"): + disp_combiner(event) + + +def test_disp_combiner_single_event_min_ang_diff(): + event_dict = { + "tel_id": [1, 2], + "hillas_intensity": [100, 200], + "hillas_width": [0.1, 0.1] * u.deg, + "hillas_length": [0.3, 0.3] * u.deg, + "hillas_fov_lon": [0.0, 0.1] * u.deg, + "hillas_fov_lat": [0.0, 0.1] * u.deg, + "hillas_psi": [5, 8] * u.deg, + "disp_tel_alt": [58.5, 58] * u.deg, + "disp_tel_az": [12.5, 15] * u.deg, + "disp_tel_parameter": [0.65, 1.1] * u.deg, + "disp_tel_sign_score": [0.95, 0.98], + "disp_tel_is_valid": [True, True], + } + event = _make_disp_event(event_dict) + + disp_combiner = StereoDispCombiner( + prefix="dummy", + min_ang_diff=10, + ) + disp_combiner(event) + + geometry = event.dl2.stereo.geometry["dummy"] + assert not geometry.is_valid + assert np.isnan(geometry.alt.to_value(u.deg)) + assert np.isnan(geometry.az.to_value(u.deg)) + + +def test_disp_combiner_n_best_tels_event(): + event_dict = { + "tel_id": [1, 2, 3], + "hillas_intensity": [150, 300, 50], + "hillas_width": [0.1, 0.2, 0.1] * u.deg, + "hillas_length": [0.3, 0.4, 0.2] * u.deg, + "hillas_fov_lon": [-0.5, 0.1, 0.3] * u.deg, + "hillas_fov_lat": [0.3, -0.2, 0.1] * u.deg, + "hillas_psi": [40, 80, -20] * u.deg, + "disp_tel_alt": [58.5, 58, 62.5] * u.deg, + "disp_tel_az": [12.5, 15, 13] * u.deg, + "disp_tel_parameter": [0.65, 1.1, 0.7] * u.deg, + "disp_tel_sign_score": [0.95, 0.98, 0.66], + "disp_tel_is_valid": [True, True, True], + } + event = _make_disp_event(event_dict) + disp_combiner = StereoDispCombiner( + prefix="dummy", + weights="intensity", + n_best_tels=2, + ) + disp_combiner(event) + + expected_event = _make_disp_event( + {key: value[:2] for key, value in event_dict.items()} + ) + expected_combiner = StereoDispCombiner( + prefix="dummy", + weights="intensity", + n_best_tels=None, + ) + expected_combiner(expected_event) + + geometry = event.dl2.stereo.geometry["dummy"] + expected_geometry = expected_event.dl2.stereo.geometry["dummy"] + assert u.isclose(geometry.alt, expected_geometry.alt) + assert u.isclose(geometry.az, expected_geometry.az) + + +def test_disp_combiner_n_tel_combinations_event(): + event_dict = { + "tel_id": [1, 2, 3], + "hillas_intensity": [150, 300, 50], + "hillas_width": [0.1, 0.2, 0.1] * u.deg, + "hillas_length": [0.3, 0.4, 0.2] * u.deg, + "hillas_fov_lon": [-0.5, 0.1, 0.3] * u.deg, + "hillas_fov_lat": [0.3, -0.2, 0.1] * u.deg, + "hillas_psi": [40, 80, -20] * u.deg, + "disp_tel_alt": [58.5, 58, 62.5] * u.deg, + "disp_tel_az": [12.5, 15, 13] * u.deg, + "disp_tel_parameter": [0.65, 1.1, 0.7] * u.deg, + "disp_tel_sign_score": [0.95, 0.98, 0.66], + "disp_tel_is_valid": [True, True, True], + } + event = _make_disp_event(event_dict) + + disp_combiner = StereoDispCombiner( + prefix="dummy", + n_tel_combinations=3, + ) + disp_combiner(event) + + stereo = event.dl2.stereo.geometry["dummy"] + assert stereo.is_valid + assert u.isclose(stereo.alt, 70.80002984 * u.deg) + assert u.isclose(stereo.az, 0.43926108 * u.deg) + + +def test_predict_disp_combiner(mono_table): + disp_combiner = StereoDispCombiner( + prefix="disp", + n_tel_combinations=2, + ) + stereo = disp_combiner.predict_table(mono_table) + + for name, field in ReconstructedGeometryContainer.fields.items(): + colname = f"disp_{name}" + assert colname in stereo.colnames + assert stereo[colname].description == field.description + + assert "obs_id" in stereo.colnames + assert "event_id" in stereo.colnames + + assert_array_equal(stereo["obs_id"], np.array([1, 1, 2, 2])) + assert_array_equal(stereo["event_id"], np.array([1, 2, 1, 2])) + assert_allclose( + stereo["disp_alt"].quantity, + [70.7338725, np.nan, 81, 70.4917615] * u.deg, + atol=1e-7, + ) + assert_allclose( + stereo["disp_az"].quantity, + [359.9419634, np.nan, 14.5, 359.5978866] * u.deg, + atol=1e-7, + ) + tel_ids = stereo["disp_telescopes"] + assert_array_equal(tel_ids[0], [1, 3]) + assert_array_equal(tel_ids[1], []) + assert_array_equal(tel_ids[2], [1]) + assert_array_equal(tel_ids[3], [1, 3, 4, 5]) + + +def test_predict_disp_combiner_n_best_tels(mono_table): + disp_combiner = StereoDispCombiner( + prefix="disp", + weights="intensity", + n_best_tels=2, + ) + stereo = disp_combiner.predict_table(mono_table) + + expected_combiner = StereoDispCombiner( + prefix="disp", + weights="intensity", + n_best_tels=None, + ) + expected = expected_combiner.predict_table(mono_table[:-2]) + + assert np.all( + u.isclose( + stereo["disp_alt"].quantity, + expected["disp_alt"].quantity, + equal_nan=True, + ) + ) + assert np.all( + u.isclose( + stereo["disp_az"].quantity, + expected["disp_az"].quantity, + equal_nan=True, + ) + ) + assert_array_equal(stereo["disp_telescopes"][-1], [1, 3]) + + +def test_predict_disp_combiner_min_ang_diff(mono_table): + disp_combiner = StereoDispCombiner( + prefix="disp", + min_ang_diff=85, + ) + stereo = disp_combiner.predict_table(mono_table) + + assert np.isnan(stereo["disp_alt"][0]) + assert np.isnan(stereo["disp_az"][0]) + assert not stereo["disp_is_valid"][0] + assert stereo["disp_telescopes"][0] == [] + + disp_combiner = StereoDispCombiner( + prefix="disp", + min_ang_diff=75, + ) + stereo = disp_combiner.predict_table(mono_table) + + assert np.isfinite(stereo["disp_alt"][0]) + assert np.isfinite(stereo["disp_az"][0]) + assert stereo["disp_is_valid"][0] + assert stereo["disp_telescopes"][0] == [1, 3] + + +def test_predict_disp_combiner_n_tel_combinations(mono_table): + disp_combiner = StereoDispCombiner( + prefix="disp", + n_tel_combinations=3, + ) + stereo = disp_combiner.predict_table(mono_table) + + assert_allclose( + stereo["disp_alt"].quantity, + [70.7338725, np.nan, 81, 70.5617748] * u.deg, + atol=1e-7, + ) + assert_allclose( + stereo["disp_az"].quantity, + [359.9419634, np.nan, 14.5, 359.84586059] * u.deg, + atol=1e-7, + ) + tel_ids = stereo["disp_telescopes"] + assert_array_equal(tel_ids[0], [1, 3]) + assert_array_equal(tel_ids[1], []) + assert_array_equal(tel_ids[2], [1]) + assert_array_equal(tel_ids[3], [1, 3, 4, 5]) + + +def test_predict_disp_combiner_empty_table(mono_table): + disp_combiner = StereoDispCombiner( + prefix="disp", + ) + empty_table = mono_table[:0] + stereo = disp_combiner.predict_table(empty_table) + cols = [ + "obs_id", + "event_id", + "disp_alt", + "disp_az", + "disp_is_valid", + "disp_telescopes", + ] + + for col in cols: + assert col in stereo.colnames + assert len(stereo) == 0 + + +def test_predict_disp_combiner_missing_disp_column(mono_table): + disp_combiner = StereoDispCombiner( + prefix="disp", + ) + broken_table = mono_table.copy() + del broken_table["disp_tel_parameter"] + with pytest.raises(KeyError, match="Required DISP column"): + disp_combiner.predict_table(broken_table) diff --git a/src/ctapipe/reco/tests/test_telescope_event_handling.py b/src/ctapipe/reco/tests/test_telescope_event_handling.py index 1faa0f4f65c..fd00163c884 100644 --- a/src/ctapipe/reco/tests/test_telescope_event_handling.py +++ b/src/ctapipe/reco/tests/test_telescope_event_handling.py @@ -1,4 +1,6 @@ +import astropy.units as u import numpy as np +import pytest from astropy.table import Table from ctapipe.io import TableLoader, read_table @@ -22,7 +24,7 @@ def test_get_subarray_index(dl1_parameters_file): ) -def test_mean_std_ufunc(dl1_parameters_file): +def test_weighted_mean_std_ufunc(dl1_parameters_file): from ctapipe.reco.telescope_event_handling import ( get_subarray_index, weighted_mean_std_ufunc, @@ -49,3 +51,259 @@ def test_mean_std_ufunc(dl1_parameters_file): assert np.allclose(mean, true_mean, equal_nan=True) assert np.allclose(std, true_std, equal_nan=True) + + +def test_get_combinations(): + from ctapipe.reco.telescope_event_handling import get_combinations + + tel_ids = [1, 2, 3] + comb_size = 2 + get_combinations.cache_clear() + index_combinations = get_combinations(len(tel_ids), comb_size) + + expected_combinations = np.array([[0, 1], [0, 2], [1, 2]]) + + assert np.allclose(index_combinations, expected_combinations) + + +def test_calc_combs_min_distances_multiple_events(): + from ctapipe.reco.telescope_event_handling import calc_combs_min_distances + + index_tel_combs = np.array( + [ + [0, 1, 2], + [3, 4, 5], + [3, 4, 6], + [3, 5, 6], + [4, 5, 6], + [7, 8, 9], + [7, 8, 10], + [7, 8, 11], + [7, 9, 10], + [7, 9, 11], + [7, 10, 11], + [8, 9, 10], + [8, 9, 11], + [8, 10, 11], + [9, 10, 11], + ] + ) + tel_ids = np.arange(12) + fov_lon_values = np.column_stack((tel_ids, tel_ids)).astype(float) + fov_lat_values = np.column_stack((tel_ids + 100, tel_ids + 100)).astype(float) + weights = tel_ids + 1 + + fov_lons, fov_lats, comb_weights = calc_combs_min_distances( + index_tel_combs, fov_lon_values, fov_lat_values, weights + ) + + expected_lons = np.array( + [ + np.average(tel_ids[combo], weights=weights[combo]) + for combo in index_tel_combs + ] + ) + expected_lats = np.array( + [ + np.average(tel_ids[combo] + 100, weights=weights[combo]) + for combo in index_tel_combs + ] + ) + expected_weights = np.array([weights[combo].sum() for combo in index_tel_combs]) + + assert np.allclose(fov_lons, expected_lons) + assert np.allclose(fov_lats, expected_lats) + assert np.allclose(comb_weights, expected_weights) + + +def test_calc_fov_lon_lat(): + from ctapipe.reco.telescope_event_handling import calc_fov_lon_lat + + prefix = "disp" + tel_table = Table( + { + "hillas_fov_lon": [1, 2, 3] * u.deg, + "hillas_fov_lat": [4, 5, 6] * u.deg, + "hillas_psi": [0, 45, 90] * u.deg, + f"{prefix}_parameter": [1, 2.5, 3] * u.deg, + } + ) + + lon, lat = calc_fov_lon_lat(tel_table, prefix) + + exp_lon = np.array([[0.0, 2.0], [0.23223305, 3.76776695], [3.0, 3.0]]) + exp_lat = np.array([[4.0, 4.0], [3.23223305, 6.76776695], [3.0, 9.0]]) + + assert np.allclose(lon, exp_lon) + assert np.allclose(lat, exp_lat) + + +def test_create_combs_array(): + from ctapipe.reco.telescope_event_handling import create_combs_array + + max_multi = 3 + k = 2 + + combs_array, combs_to_multi_indices = create_combs_array(max_multi, k) + + exp_combs_array = np.array([[0, 1], [0, 1], [0, 2], [1, 2]]) + exp_combs_to_multi_indices = np.array([2, 3, 3, 3]) + + assert np.allclose(combs_array, exp_combs_array) + assert np.allclose(combs_to_multi_indices, exp_combs_to_multi_indices) + + +def test_get_index_combs(): + from ctapipe.reco.telescope_event_handling import get_index_combs + + multiplicities = np.array([2, 1, 3, 2]) + combs_array = np.array([[0, 1], [0, 1], [0, 2], [1, 2]]) + combs_to_multi_indices = np.array([2, 3, 3, 3]) + k = 2 + + index_tel_combs, num_combs = get_index_combs( + multiplicities, combs_array, combs_to_multi_indices, k + ) + + exp_index_tel_combs = np.array([[0, 1], [3, 4], [3, 5], [4, 5], [6, 7]]) + exp_num_combs = np.array([1, 0, 3, 1]) + + assert np.allclose(index_tel_combs, exp_index_tel_combs) + assert np.allclose(num_combs, exp_num_combs) + + +@pytest.mark.parametrize( + ("k", "expected"), + [ + (2, np.array([[0, 0], [0, 1], [1, 0], [1, 1]])), + (0, np.empty((1, 0), dtype=int)), + ], +) +def test_binary_combinations(k, expected): + from ctapipe.reco.telescope_event_handling import binary_combinations + + combinations = binary_combinations(k) + + assert combinations.shape == expected.shape + assert np.array_equal(combinations, expected) + + +@pytest.mark.parametrize( + ("min_ang_diff", "psi1", "psi2", "expected"), + [ + (20.0, 10.0 * u.deg, 50.0 * u.deg, True), + (10.0, 0.0 * u.deg, 180.0 * u.deg, False), + ], +) +def test_check_ang_diff(min_ang_diff, psi1, psi2, expected): + from ctapipe.reco.telescope_event_handling import check_ang_diff + + result = check_ang_diff(min_ang_diff, psi1, psi2) + + assert bool(result) is expected + + +def test_check_ang_diff_array(): + from ctapipe.reco.telescope_event_handling import check_ang_diff + + min_ang_diff = 20.0 + psi1 = np.array([0.0, 10.0, 0.0]) * u.deg + psi2 = np.array([30.0, 25.0, 180.0]) * u.deg + + result = check_ang_diff(min_ang_diff, psi1, psi2) + expected = np.array([True, False, False]) + + assert np.all(result == expected) + + +def test_valid_tels_of_multi(): + from ctapipe.reco.telescope_event_handling import valid_tels_of_multi + + valid_tel_to_array_indices = np.array([0, 0, 1, 1, 1, 2]) + + mask = valid_tels_of_multi(2, valid_tel_to_array_indices) + expected = np.array([True, True, False, False, False, False]) + + assert np.array_equal(mask, expected) + + +def test_valid_tels_of_multi_noop(): + from ctapipe.reco.telescope_event_handling import valid_tels_of_multi + + valid_tel_to_array_indices = np.array([0, 0, 1, 1, 1, 2]) + + mask = valid_tels_of_multi(4, valid_tel_to_array_indices) + + assert not np.any(mask) + + +def test_fill_lower_multiplicities(): + from ctapipe.reco.telescope_event_handling import fill_lower_multiplicities + + fov_lon_combs_mean = np.array([10.0, np.nan, 30.0]) + fov_lat_combs_mean = np.array([0.0, np.nan, 1.0]) + n_tel_combinations = 3 + valid_tel_to_array_indices = np.array([0, 0, 0, 1, 1, 2]) + valid_multiplicity = np.array([3, 2, 1]) + fov_lon_values = np.array( + [ + [0.0, 1.0], + [0.0, 1.0], + [0.0, 1.0], + [1.0, 11.0], + [2.0, 12.0], + [5.0, 15.0], + ] + ) + fov_lat_values = np.array( + [ + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [1.0, 1.0], + ] + ) + weights = np.ones(len(valid_tel_to_array_indices)) + + fill_lower_multiplicities( + fov_lon_combs_mean, + fov_lat_combs_mean, + n_tel_combinations, + valid_tel_to_array_indices, + valid_multiplicity, + fov_lon_values, + fov_lat_values, + weights, + ) + + assert np.allclose(fov_lon_combs_mean, np.array([10.0, 1.5, 30.0]), equal_nan=True) + assert np.allclose(fov_lat_combs_mean, np.array([0.0, 0.0, 1.0]), equal_nan=True) + + +def test_fill_lower_multiplicities_noop(): + from ctapipe.reco.telescope_event_handling import fill_lower_multiplicities + + fov_lon_combs_mean = np.array([10.0, 20.0]) + fov_lat_combs_mean = np.array([0.0, 1.0]) + n_tel_combinations = 2 + valid_tel_to_array_indices = np.array([0, 0, 1, 1]) + valid_multiplicity = np.array([2, 2]) + fov_lon_values = np.array([[0.0, 1.0], [0.0, 1.0], [2.0, 3.0], [2.0, 3.0]]) + fov_lat_values = np.zeros_like(fov_lon_values) + weights = np.ones(len(valid_tel_to_array_indices)) + + fill_lower_multiplicities( + fov_lon_combs_mean, + fov_lat_combs_mean, + n_tel_combinations, + valid_tel_to_array_indices, + valid_multiplicity, + fov_lon_values, + fov_lat_values, + weights, + ) + + assert np.allclose(fov_lon_combs_mean, np.array([10.0, 20.0])) + assert np.allclose(fov_lat_combs_mean, np.array([0.0, 1.0])) diff --git a/src/ctapipe/reco/utils.py b/src/ctapipe/reco/utils.py index e4327094dcd..79fcb23106f 100644 --- a/src/ctapipe/reco/utils.py +++ b/src/ctapipe/reco/utils.py @@ -1,3 +1,17 @@ +import astropy.units as u +import numpy as np + + +def _default_column_for_length(default, n_rows): + """Create a length-matching column filled with `default` (unit-safe).""" + if n_rows == 0: + if isinstance(default, u.Quantity): + return u.Quantity([], unit=default.unit) + return np.full(0, default) + + return default + + def add_defaults_and_meta(table, container, prefix=None, add_tel_prefix=False): """ Fill column descriptions and default values into table for container @@ -17,17 +31,17 @@ def add_defaults_and_meta(table, container, prefix=None, add_tel_prefix=False): if prefix is None: prefix = container.default_prefix + n_rows = len(table) + col_prefix = f"{prefix}_tel_" if add_tel_prefix else f"{prefix}_" + for name, field in container.fields.items(): if add_tel_prefix and name == "telescopes": continue - if add_tel_prefix: - colname = f"{prefix}_tel_{name}" - else: - colname = f"{prefix}_{name}" + colname = f"{col_prefix}{name}" if colname not in table.colnames and field.default is not None: - table[colname] = field.default + table[colname] = _default_column_for_length(field.default, n_rows) if colname in table.colnames: table[colname].description = field.description diff --git a/src/ctapipe/resources/train_disp_reconstructor.yaml b/src/ctapipe/resources/train_disp_reconstructor.yaml index 47a43ead025..ba6acadd916 100644 --- a/src/ctapipe/resources/train_disp_reconstructor.yaml +++ b/src/ctapipe/resources/train_disp_reconstructor.yaml @@ -6,8 +6,8 @@ # ====================================================================== TrainDispReconstructor: - random_seed: 0 # Seed used for sampling n_events for training. - n_events: # The number of events used for training that can be provided + random_seed: 0 # Seed used for sampling n_events for training. + n_events: # The number of events used for training that can be provided # - [type, "LST*", 1000] # independently for each telescope type (e.g. "LST_LST_LSTCam"). # - [type, "MST*", 1000] # If not specified, all events in the file are used. @@ -17,7 +17,7 @@ TrainDispReconstructor: DispReconstructor: # prefix: # Add a prefix of the output here, if you want to apply multiple - # DispReconstructors on the same file (e.g. for comparing different settings) + # DispReconstructors on the same file (e.g. for comparing different settings) # How many cores to use. Overwrites model config n_jobs: -1 @@ -38,18 +38,18 @@ TrainDispReconstructor: n_estimators: 10 max_depth: 10 - QualityQuery: # Event Selection performed before training the models + MLQualityQuery: # Event Selection performed before training the models quality_criteria: - # - [, ] + # - [, ] - ["HillasValid", "HillasReconstructor_is_valid"] - ["enough intensity", "hillas_intensity > 50"] - ["Positive width", "hillas_width > 0"] - ["enough pixels", "morphology_n_pixels > 3"] - ["not clipped", "leakage_intensity_width_2 < 0.5"] - FeatureGenerator: # On-the-fly generation of additional features + FeatureGenerator: # On-the-fly generation of additional features features: - # - [, ] + # - [, ] - ["area", "hillas_width * hillas_length"] features: @@ -80,10 +80,10 @@ TrainDispReconstructor: - intensity_std - intensity_skewness - intensity_kurtosis - - ExtraTreesRegressor_energy # Mean predicted energy of the array event - - ExtraTreesRegressor_tel_energy # Energy prediction for single telescope image - - area # Features generated by the FeatureGenerator still have to be listed here - # for the models to use them. + - ExtraTreesRegressor_energy # Mean predicted energy of the array event + - ExtraTreesRegressor_tel_energy # Energy prediction for single telescope image + - area # Features generated by the FeatureGenerator still have to be listed here + # for the models to use them. # stereo_combiner_cls: "StereoMeanCombiner" # StereoMeanCombiner: # Here you can set options for StereoMeanCombiner diff --git a/src/ctapipe/resources/train_energy_regressor.yaml b/src/ctapipe/resources/train_energy_regressor.yaml index 61778c3181b..8149ce6d9e5 100644 --- a/src/ctapipe/resources/train_energy_regressor.yaml +++ b/src/ctapipe/resources/train_energy_regressor.yaml @@ -6,8 +6,8 @@ # ====================================================================== TrainEnergyRegressor: - random_seed: 0 # Seed used for sampling n_events for training. - n_events: # The number of events used for training that can be provided + random_seed: 0 # Seed used for sampling n_events for training. + n_events: # The number of events used for training that can be provided # - [type, "LST*", 1000] # independently for each telescope type (e.g. "LST_LST_LSTCam"). # - [type, "MST*", 1000] # If not specified, all events in the file are used. @@ -17,7 +17,7 @@ TrainEnergyRegressor: EnergyRegressor: # prefix: # Add a prefix of the output here, if you want to apply multiple - # EnergyRegressors on the same file (e.g. for comparing different settings) + # EnergyRegressors on the same file (e.g. for comparing different settings) log_target: True n_jobs: -1 @@ -29,21 +29,24 @@ TrainEnergyRegressor: n_estimators: 10 max_depth: 10 - QualityQuery: # Event Selection performed before training the models + MLQualityQuery: # Event Selection performed before training the models quality_criteria: - # - [, ] + # - [, ] - ["HillasValid", "HillasReconstructor_is_valid"] - ["enough intensity", "hillas_intensity > 50"] - ["Positive width", "hillas_width > 0"] - ["enough pixels", "morphology_n_pixels > 3"] - ["not clipped", "leakage_intensity_width_2 < 0.5"] - FeatureGenerator: # On-the-fly generation of additional features + FeatureGenerator: # On-the-fly generation of additional features features: - # - [, ] + # - [, ] - ["area", "hillas_width * hillas_length"] - ["n_telescopes_triggered", "subarray.multiplicity(tels_with_trigger)"] - - ["n_telescopes_hillas_reconstructor", "subarray.multiplicity(HillasReconstructor_telescopes)"] + - [ + "n_telescopes_hillas_reconstructor", + "subarray.multiplicity(HillasReconstructor_telescopes)", + ] features: - hillas_intensity @@ -76,9 +79,10 @@ TrainEnergyRegressor: - intensity_std - intensity_skewness - intensity_kurtosis - - area # Features generated by the FeatureGenerator - - n_telescopes_triggered # still have to be listed here for the models - - n_telescopes_hillas_reconstructor # to use them. + - area # Features generated by the FeatureGenerator + - n_telescopes_triggered # still have to be listed here for the models + - n_telescopes_hillas_reconstructor # to use them. + # stereo_combiner_cls: "StereoMeanCombiner" # StereoMeanCombiner: # Here you can set options for StereoMeanCombiner diff --git a/src/ctapipe/resources/train_particle_classifier.yaml b/src/ctapipe/resources/train_particle_classifier.yaml index 5a33d72945e..1b812eaa1bf 100644 --- a/src/ctapipe/resources/train_particle_classifier.yaml +++ b/src/ctapipe/resources/train_particle_classifier.yaml @@ -6,11 +6,11 @@ # ====================================================================== TrainParticleClassifier: - random_seed: 0 # Seed used for sampling n_events for training. - n_events: # The number of events used for training that can be provided + random_seed: 0 # Seed used for sampling n_events for training. + n_events: # The number of events used for training that can be provided # - [type, "LST*", 1000] # independently for each telescope type (e.g. "LST_LST_LSTCam"). # - [type, "MST*", 1000] # If not specified, as many events as possible are used. - signal_fraction: 0.5 # signal_fraction = n_signal / n_events + signal_fraction: 0.5 # signal_fraction = n_signal / n_events CrossValidator: n_cross_validations: 5 @@ -18,7 +18,7 @@ TrainParticleClassifier: ParticleClassifier: # prefix: # Add a prefix of the output here, if you want to apply multiple - # ParticleClassifiers on the same file (e.g. for comparing different settings) + # ParticleClassifiers on the same file (e.g. for comparing different settings) # How many cores to use. Overwrites model config n_jobs: -1 @@ -30,18 +30,18 @@ TrainParticleClassifier: n_estimators: 10 max_depth: 10 - QualityQuery: # Event Selection performed before training the models + MLQualityQuery: # Event Selection performed before training the models quality_criteria: - # - [, ] + # - [, ] - ["HillasValid", "HillasReconstructor_is_valid"] - ["enough intensity", "hillas_intensity > 50"] - ["Positive width", "hillas_width > 0"] - ["enough pixels", "morphology_n_pixels > 3"] - ["not clipped", "leakage_intensity_width_2 < 0.5"] - FeatureGenerator: # On-the-fly generation of additional features + FeatureGenerator: # On-the-fly generation of additional features features: - # - [, ] + # - [, ] - ["area", "hillas_width * hillas_length"] features: @@ -75,10 +75,10 @@ TrainParticleClassifier: - intensity_std - intensity_skewness - intensity_kurtosis - - ExtraTreesRegressor_energy # Mean predicted energy of the array event - - ExtraTreesRegressor_tel_energy # Energy prediction for single telescope image - - area # Features generated by the FeatureGenerator still have to be listed here - # for the models to use them. + - ExtraTreesRegressor_energy # Mean predicted energy of the array event + - ExtraTreesRegressor_tel_energy # Energy prediction for single telescope image + - area # Features generated by the FeatureGenerator still have to be listed here + # for the models to use them. # stereo_combiner_cls: "StereoMeanCombiner" # StereoMeanCombiner: # Here you can set options for StereoMeanCombiner diff --git a/src/ctapipe/tests/test_traitlets_configurable.py b/src/ctapipe/tests/test_traitlets_configurable.py index 696e3dac831..78b2ed24347 100644 --- a/src/ctapipe/tests/test_traitlets_configurable.py +++ b/src/ctapipe/tests/test_traitlets_configurable.py @@ -29,6 +29,10 @@ "classes", } +not_configurable_traits = { + "ctapipe.reco.stereo_combination.StereoDispCombiner": {"property"}, +} + def test_all_traitlets_configurable(): skip_modules = {"_dev_version", "tests"} @@ -64,14 +68,15 @@ def find_all_traitlets(module, missing_config=None): has_traits = False if has_traits: + fqcn = f"{obj.__module__}.{obj.__name__}" + exempt = not_configurable_traits.get(fqcn, set()) + for traitname, trait in obj.class_traits().items(): - if ( - not trait.metadata.get("config", False) - and traitname not in ignore_traits - ): - missing_config[f"{obj.__module__}.{obj.__name__}"].add( - traitname - ) + if traitname in exempt or traitname in ignore_traits: + continue + + if not trait.metadata.get("config", False): + missing_config[fqcn].add(traitname) return missing_config