From f3d441d01465c8e2e369ee664b9b5857a964eaf9 Mon Sep 17 00:00:00 2001 From: gschwefer Date: Tue, 25 Feb 2025 12:07:08 +0100 Subject: [PATCH 01/43] Include template paths as telescope parameters --- src/ctapipe/reco/impact.py | 203 +++++++----------- .../utils/template_network_interpolator.py | 40 +++- 2 files changed, 114 insertions(+), 129 deletions(-) diff --git a/src/ctapipe/reco/impact.py b/src/ctapipe/reco/impact.py index 9f021d3c56a..0e3ba702eff 100644 --- a/src/ctapipe/reco/impact.py +++ b/src/ctapipe/reco/impact.py @@ -3,7 +3,6 @@ Implementation of the ImPACT reconstruction algorithm """ import copy -from string import Template import numpy as np import numpy.ma as ma @@ -12,6 +11,7 @@ from scipy.stats import norm from ctapipe.core import traits +from ctapipe.core.telescope_component import TelescopeParameter from ctapipe.exceptions import OptionalDependencyMissing from ..compat import COPY_IF_NEEDED @@ -31,8 +31,6 @@ neg_log_likelihood_approx, ) from ..utils.template_network_interpolator import ( - DummyTemplateInterpolator, - DummyTimeInterpolator, TemplateNetworkInterpolator, TimeGradientInterpolator, ) @@ -109,13 +107,18 @@ class ImPACTReconstructor(HillasGeometryReconstructor): """ - use_time_gradient = traits.Bool( - default_value=False, - help="Use time gradient in ImPACT reconstruction. Requires an extra set of time gradient templates", + image_template_path = TelescopeParameter( + trait=traits.Path(exists=True, directory_ok=False, allow_none=False), + allow_none=False, + help=("Path to the image templates to be used in the reconstruction"), ).tag(config=True) - root_dir = traits.Unicode( - default_value=".", help="Directory containing ImPACT tables" + # The time gradient templates are optional, so None is allowed here + time_gradient_template_path = TelescopeParameter( + trait=traits.Path(exists=True, directory_ok=False, allow_none=True), + allow_none=True, + default_value=None, + help=("Path to the time gradient templates to be used in the reconstruction"), ).tag(config=True) # For likelihood calculation we need the with of the @@ -149,14 +152,12 @@ def __init__( super().__init__(subarray, atmosphere_profile, **kwargs) + self.subarray = subarray + # First we create a dictionary of image template interpolators # for each telescope type # self.priors = prior - # String templates for loading ImPACT templates - self.amplitude_template = Template("${base}/${camera}.template.gz") - self.time_template = Template("${base}/${camera}_time.template.gz") - # Next we need the position, area and amplitude from each pixel in the event # making this a class member makes passing them around much easier @@ -263,53 +264,61 @@ def __call__(self, event): self._store_impact_parameter(event) - def initialise_templates(self, tel_type): - """Check if templates for a given telescope type has been initialised - and if not do it and add to the dictionary + def set_up_templates(self, tel_ids): + template_sort_dict = {} + time_template_sort_dict = {} - Parameters - ---------- - tel_type: dictionary - Dictionary of telescope types in event - - Returns - ------- - boolean: Confirm initialisation + self.use_time_gradient = True - """ + for tel_id in tel_ids: + if self.image_template_path.tel[tel_id] not in template_sort_dict.keys(): + template_sort_dict[self.image_template_path.tel[tel_id]] = [tel_id] + else: + template_sort_dict[self.image_template_path.tel[tel_id]].append(tel_id) - for t in tel_type: - if tel_type[t] in self.prediction.keys() or tel_type[t] == "DUMMY": - continue + if self.time_gradient_template_path.tel[tel_id] is not None: + if ( + self.time_gradient_template_path.tel[tel_id] + not in time_template_sort_dict.keys() + ): + time_template_sort_dict[ + self.time_gradient_template_path.tel[tel_id] + ] = [tel_id] + else: + time_template_sort_dict[ + self.time_gradient_template_path.tel[tel_id] + ].append(tel_id) - if self.dummy_reconstructor: - self.prediction[tel_type[t]] = DummyTemplateInterpolator() else: - filename = self.amplitude_template.substitute( - base=self.root_dir, camera=tel_type[t] - ) - self.prediction[tel_type[t]] = TemplateNetworkInterpolator( - filename, bounds=((-5, 1), (-1.5, 1.5)) - ) - PROV.add_input_file( - filename, role="ImPACT Template file for " + tel_type[t] - ) + self.use_time_gradient = False - if self.use_time_gradient: - if self.dummy_reconstructor: - self.time_prediction[tel_type[t]] = DummyTimeInterpolator() - else: - filename = self.time_template.substitute( - base=self.root_dir, camera=tel_type[t] - ) - self.time_prediction[tel_type[t]] = TimeGradientInterpolator( - filename - ) - PROV.add_input_file( - filename, role="ImPACT Time Template file for " + tel_type[t] + for template_path, tel_ids in template_sort_dict.items(): + net_interpolator = TemplateNetworkInterpolator( + template_path, bounds=((-5, 1), (-1.5, 1.5)) + ) + + interp_tel_string = net_interpolator.tel_type_string + for id in tel_ids: + if interp_tel_string != str(self.subarray.tel[id]): + raise ValueError( + "You are using templates that are not intended for this telescope type" ) - return True + self.prediction[tuple(tel_ids)] = net_interpolator + + if self.use_time_gradient: + for template_path, tel_ids in time_template_sort_dict.items(): + time_interpolator = TimeGradientInterpolator(template_path) + + interp_tel_string = time_interpolator.tel_type_string + + for id in tel_ids: + if interp_tel_string != str(self.subarray.tel[id]): + raise ValueError( + "You are using templates that are not intended for this telescope type" + ) + + self.time_prediction[tuple(tel_ids)] = time_interpolator def get_hillas_mean(self): """This is a simple function to find the peak position of each image @@ -400,59 +409,6 @@ def get_shower_max(self, source_x, source_y, core_x, core_y, zen): return slant_depth.to_value(u.g / (u.cm * u.cm)) - def image_prediction( - self, tel_type, zenith, azimuth, energy, impact, x_max, pix_x, pix_y - ): - """Creates predicted image for the specified pixels, interpolated - from the template library. - - Parameters - ---------- - tel_type: string - Telescope type specifier - energy: float - Event energy (TeV) - impact: float - Impact diance of shower (metres) - x_max: float - Depth of shower maximum (num bins from expectation) - pix_x: ndarray - X coordinate of pixels - pix_y: ndarray - Y coordinate of pixels - - Returns - ------- - ndarray: predicted amplitude for all pixels - - """ - return self.prediction[tel_type]( - zenith, azimuth, energy, impact, x_max, pix_x, pix_y - ) - - def predict_time(self, tel_type, zenith, azimuth, energy, impact, x_max): - """Creates predicted image for the specified pixels, interpolated - from the template library. - - Parameters - ---------- - tel_type: string - Telescope type specifier - energy: float - Event energy (TeV) - impact: float - Impact diance of shower (metres) - x_max: float - Depth of shower maximum (num bins from expectation) - - Returns - ------- - ndarray: predicted amplitude for all pixels - - """ - time = self.time_prediction[tel_type](zenith, azimuth, energy, impact, x_max) - return time.T[0], time.T[1] - def get_likelihood( self, source_x, @@ -540,33 +496,34 @@ def get_likelihood( np.zeros(self.image.shape[0]), ) # Loop over all telescope types and get prediction - for tel_type in np.unique(self.tel_types).tolist(): - type_mask = self.tel_types == tel_type - prediction[type_mask] = self.image_prediction( - tel_type, + for tel_ids, template in self.prediction.items(): + template_mask = np.array([id in tel_ids for id in self.tel_id]) + + prediction[template_mask] = template( np.rad2deg(zenith), azimuth, - energy * np.ones_like(impact[type_mask]), - impact[type_mask], - x_max_diff * np.ones_like(impact[type_mask]), - np.rad2deg(pix_x_rot[type_mask]), - np.rad2deg(pix_y_rot[type_mask]), + energy * np.ones_like(impact[template_mask]), + impact[template_mask], + x_max_diff * np.ones_like(impact[template_mask]), + np.rad2deg(pix_x_rot[template_mask]), + np.rad2deg(pix_y_rot[template_mask]), ) - if self.use_time_gradient: - tg, tgu = self.predict_time( - tel_type, + if self.use_time_gradient: + for tel_ids, time_template in self.time_prediction.items(): + time_template_mask = np.array([id in tel_ids for id in self.tel_id]) + time_pred = time_template( np.rad2deg(zenith), azimuth, - energy * np.ones_like(impact[type_mask]), - impact[type_mask], - x_max_diff * np.ones_like(impact[type_mask]), + energy * np.ones_like(impact[time_template_mask]), + impact[time_template_mask], + x_max_diff * np.ones_like(impact[time_template_mask]), ) - time_gradients[type_mask] = tg - time_gradients_uncertainty[type_mask] = tgu - if self.use_time_gradient: + time_gradients[time_template_mask] = time_pred.T[0] + time_gradients_uncertainty[time_template_mask] = time_pred.T[1] + time_gradients_uncertainty[time_gradients_uncertainty == 0] = 1e-6 chi2 = 0 @@ -765,7 +722,7 @@ class before minimisation can take place. This simply copies a # Finally run some functions to get ready for the event self.get_hillas_mean() - self.initialise_templates(type_tel) + self.set_up_templates(hillas_dict.keys()) def predict( self, diff --git a/src/ctapipe/utils/template_network_interpolator.py b/src/ctapipe/utils/template_network_interpolator.py index e921056d280..78578af06c5 100644 --- a/src/ctapipe/utils/template_network_interpolator.py +++ b/src/ctapipe/utils/template_network_interpolator.py @@ -269,12 +269,21 @@ def __init__(self, template_file, bounds=((-5, 1), (-1.5, 1.5))): super().__init__() + self.template_file = template_file + input_dict = None with gzip.open(template_file, "r") as file_list: input_dict = pickle.load(file_list) - keys = np.array(list(input_dict.keys())) - values = np.array(list(input_dict.values()), dtype=np.float32) + # The dict has at its highest level two keys: data and tel_type. + # The first just contains the template dict, the second a string to validate the tel_type + # that the templates are made for. + + data_input_dict = input_dict["data"] + self.tel_type_string = input_dict["tel_type"] + + keys = np.array(list(data_input_dict.keys())) + values = np.array(list(data_input_dict.values()), dtype=np.float32) self.no_zenaz = False self.bounds = bounds @@ -324,6 +333,10 @@ def __call__(self, zenith, azimuth, energy, impact, xmax, xb, yb): return interpolated_value + # Two interpolators are the same if they are generated from the same file + def __eq__(self, other): + return self.template_file == other.template_file + class TimeGradientInterpolator(BaseTemplate): """ @@ -339,11 +352,22 @@ def __init__(self, template_file): """ super().__init__() - file_list = gzip.open(template_file) - input_dict = pickle.load(file_list) - keys = np.array(list(input_dict.keys())) - values = np.array(list(input_dict.values()), dtype=np.float32) + self.template_file = template_file + + input_dict = None + with gzip.open(template_file, "r") as file_list: + input_dict = pickle.load(file_list) + + # The dict has at its highest level two keys: data and tel_type. + # The first just contains the template dict, the second a string to validate the tel_type + # that the templates are made for. + + data_input_dict = input_dict["data"] + self.tel_type_string = input_dict["tel_type"] + + keys = np.array(list(data_input_dict.keys())) + values = np.array(list(data_input_dict.values()), dtype=np.float32) self.no_zenaz = False # First check if we even have a zen and azimuth entry @@ -385,6 +409,10 @@ def __call__(self, zenith, azimuth, energy, impact, xmax): return interpolated_value + # Two interpolators are the same if they are generated from the same file + def __eq__(self, other): + return self.template_file == other.template_file + class DummyTemplateInterpolator: def __call__(self, zenith, azimuth, energy, impact, xmax, xb, yb): From cd6e775d8b3d63b33d447c9e177b212227165058 Mon Sep 17 00:00:00 2001 From: gschwefer Date: Tue, 25 Feb 2025 18:23:57 +0100 Subject: [PATCH 02/43] Move template set up to reconstructor init --- src/ctapipe/reco/impact.py | 59 ++++++++++++++++++++++---------------- 1 file changed, 35 insertions(+), 24 deletions(-) diff --git a/src/ctapipe/reco/impact.py b/src/ctapipe/reco/impact.py index 0e3ba702eff..d2996825647 100644 --- a/src/ctapipe/reco/impact.py +++ b/src/ctapipe/reco/impact.py @@ -176,6 +176,8 @@ def __init__( self.prediction = dict() self.time_prediction = dict() + self.set_up_templates() + self.array_direction = None self.nominal_frame = None @@ -264,13 +266,13 @@ def __call__(self, event): self._store_impact_parameter(event) - def set_up_templates(self, tel_ids): + def set_up_templates(self): template_sort_dict = {} time_template_sort_dict = {} self.use_time_gradient = True - for tel_id in tel_ids: + for tel_id in self.subarray.tel_ids: if self.image_template_path.tel[tel_id] not in template_sort_dict.keys(): template_sort_dict[self.image_template_path.tel[tel_id]] = [tel_id] else: @@ -498,31 +500,32 @@ def get_likelihood( # Loop over all telescope types and get prediction for tel_ids, template in self.prediction.items(): - template_mask = np.array([id in tel_ids for id in self.tel_id]) - - prediction[template_mask] = template( - np.rad2deg(zenith), - azimuth, - energy * np.ones_like(impact[template_mask]), - impact[template_mask], - x_max_diff * np.ones_like(impact[template_mask]), - np.rad2deg(pix_x_rot[template_mask]), - np.rad2deg(pix_y_rot[template_mask]), - ) - - if self.use_time_gradient: - for tel_ids, time_template in self.time_prediction.items(): - time_template_mask = np.array([id in tel_ids for id in self.tel_id]) - time_pred = time_template( + template_mask = self.template_masks[tel_ids] + if np.any(template_mask): + prediction[template_mask] = template( np.rad2deg(zenith), azimuth, - energy * np.ones_like(impact[time_template_mask]), - impact[time_template_mask], - x_max_diff * np.ones_like(impact[time_template_mask]), + energy * np.ones_like(impact[template_mask]), + impact[template_mask], + x_max_diff * np.ones_like(impact[template_mask]), + np.rad2deg(pix_x_rot[template_mask]), + np.rad2deg(pix_y_rot[template_mask]), ) - time_gradients[time_template_mask] = time_pred.T[0] - time_gradients_uncertainty[time_template_mask] = time_pred.T[1] + if self.use_time_gradient: + for tel_ids, time_template in self.time_prediction.items(): + time_template_mask = self.time_template_masks[tel_ids] + if np.any(time_template_mask): + time_pred = time_template( + np.rad2deg(zenith), + azimuth, + energy * np.ones_like(impact[time_template_mask]), + impact[time_template_mask], + x_max_diff * np.ones_like(impact[time_template_mask]), + ) + + time_gradients[time_template_mask] = time_pred.T[0] + time_gradients_uncertainty[time_template_mask] = time_pred.T[1] time_gradients_uncertainty[time_gradients_uncertainty == 0] = 1e-6 @@ -720,9 +723,17 @@ class before minimisation can take place. This simply copies a self.image[mask] = ma.masked self.time[mask] = ma.masked + self.template_masks = { + tels_key: np.isin(list(hillas_dict.keys()), tels_key) + for tels_key in self.prediction.keys() + } + if self.use_time_gradient: + self.time_template_masks = { + tels_key: np.isin(list(hillas_dict.keys()), tels_key) + for tels_key in self.time_prediction.keys() + } # Finally run some functions to get ready for the event self.get_hillas_mean() - self.set_up_templates(hillas_dict.keys()) def predict( self, From 5daf0cd89375fd279925fcee6fc0f0e3c5f598dc Mon Sep 17 00:00:00 2001 From: gschwefer Date: Fri, 28 Feb 2025 14:00:20 +0100 Subject: [PATCH 03/43] Make spe and pedestal width telescope parameters --- src/ctapipe/reco/impact.py | 100 +++++++++++++++++++++++++++---------- 1 file changed, 75 insertions(+), 25 deletions(-) diff --git a/src/ctapipe/reco/impact.py b/src/ctapipe/reco/impact.py index d2996825647..13200d699bb 100644 --- a/src/ctapipe/reco/impact.py +++ b/src/ctapipe/reco/impact.py @@ -69,6 +69,29 @@ MINUIT_STRATEGY = 1 # Default minimization strategy, 2 is careful, 0 is fast MINUIT_TOLERANCE_FACTOR = 1000 # Tolerance for convergence according to EDM criterion MIGRAD_ITERATE = 1 # Do not call migrad again if convergence was not reached + +BACKUP_PED_TABLE = { + "LST_LST_LSTCam": 1.4, + "MST_MST_NectarCam": 1.3, + "MST_MST_FlashCam": 1.3, + "SST_SST_SST-Camera": 0.5, + "SST_SST_CHEC": 0.5, + "SST_ASTRI_ASTRICam": 0.5, + "dummy": 0.01, + "UNKNOWN-960PX": 1.0, +} + +BACKUP_SPE_TABLE = { + "LST_LST_LSTCam": 0.6, + "MST_MST_NectarCam": 0.6, + "MST_MST_FlashCam": 0.6, + "SST_SST_SST-Camera": 0.6, + "SST_SST_CHEC": 0.6, + "SST_ASTRI_ASTRICam": 0.6, + "dummy": 0.6, + "UNKNOWN-960PX": 0.6, +} + __all__ = ["ImPACTReconstructor"] @@ -121,21 +144,21 @@ class ImPACTReconstructor(HillasGeometryReconstructor): help=("Path to the time gradient templates to be used in the reconstruction"), ).tag(config=True) - # For likelihood calculation we need the with of the - # pedestal distribution for each pixel - # currently this is not available from the calibration, - # so for now lets hard code it in a dict - ped_table = { - "LSTCam": 1.4, - "NectarCam": 1.3, - "FlashCam": 1.3, - "SST-Camera": 0.5, - "CHEC": 0.5, - "ASTRICam": 0.5, - "dummy": 0.01, - "UNKNOWN-960PX": 1.0, - } - spe = 0.6 # Also hard code single p.e. distribution width + # The SPE and pedestal width parameters are also configurable as TelescopeParameters. + # None is allowed and the default value. In that case, either values from the event monitoring data are used + # or if that is also not available, a value from a hardcoded backup dict is used. + + pedestal_width = traits.FloatTelescopeParameter( + allow_none=True, + default_value=None, + help="Pedestal width parameter for the likelihood", + ).tag(config=True) + + spe_width = traits.FloatTelescopeParameter( + allow_none=True, + default_value=None, + help="SPE width parameter for the likelihood", + ).tag(config=True) property = ReconstructionProperty.ENERGY | ReconstructionProperty.GEOMETRY @@ -164,7 +187,7 @@ def __init__( self.pixel_x, self.pixel_y = None, None self.image, self.time = None, None - self.tel_types, self.tel_id = None, None + self.tel_id = None # We also need telescope positions self.tel_pos_x, self.tel_pos_y = None, None @@ -212,12 +235,32 @@ def __call__(self, event): telescope_pointings = self._get_telescope_pointings(event) # Finally get the telescope images and and the selection masks - mask_dict, image_dict, time_dict = {}, {}, {} + mask_dict, image_dict, time_dict, ped_dict, spe_dict = {}, {}, {}, {}, {} for tel_id in hillas_dict.keys(): image = event.dl1.tel[tel_id].image image_dict[tel_id] = image time_dict[tel_id] = event.dl1.tel[tel_id].peak_time mask = event.dl1.tel[tel_id].image_mask + if self.pedestal_width.tel[tel_id] is None: + if event.mon.tel[tel_id].pedestal.charge_std is None: + ped_dict[tel_id] = BACKUP_PED_TABLE[ + str(self.subarray.tel[tel_id]) + ] * np.ones(self.subarray.tel[tel_id].camera.geometry.n_pixels) + else: + ped_dict[tel_id] = event.mon.tel[tel_id].pedestal.charge_std[0] + else: + ped_dict[tel_id] = self.pedestal_width.tel[tel_id] * np.ones( + self.subarray.tel[tel_id].camera.geometry.n_pixels + ) + + if self.spe_width.tel[tel_id] is None: + spe_dict[tel_id] = BACKUP_SPE_TABLE[ + str(self.subarray.tel[tel_id]) + ] * np.ones(self.subarray.tel[tel_id].camera.geometry.n_pixels) + else: + spe_dict[tel_id] = self.spe_width.tel[tel_id] * np.ones( + self.subarray.tel[tel_id].camera.geometry.n_pixels + ) # Dilate the images around the original cleaning to help the fit for _ in range(3): @@ -260,6 +303,8 @@ def __call__(self, event): image_dict=image_dict, mask_dict=mask_dict, time_dict=time_dict, + ped_dict=ped_dict, + spe_dict=spe_dict, ) event.dl2.stereo.geometry[self.__class__.__name__] = shower_result event.dl2.stereo.energy[self.__class__.__name__] = energy_result @@ -591,6 +636,8 @@ def set_event_properties( image_dict, time_dict, mask_dict, + ped_dict, + spe_dict, subarray, array_pointing, telescope_pointing, @@ -630,11 +677,10 @@ class before minimisation can take place. This simply copies a self.tel_pos_y = np.zeros(len(hillas_dict)) # self.scale_factor = np.zeros(len(hillas_dict)) - self.ped = np.zeros(len(hillas_dict)) - self.tel_types, self.tel_id = list(), list() + self.tel_id = list() max_pix_x = 0 - px, py, pa, pt = list(), list(), list(), list() + px, py, pa, pt, p_ped, p_spe = list(), list(), list(), list(), list(), list() self.hillas_parameters = list() # Get telescope positions in tilted frame @@ -682,8 +728,9 @@ class before minimisation can take place. This simply copies a pa.append(image_dict[tel_id][mask]) pt.append(time_dict[tel_id][mask]) - self.ped[i] = self.ped_table[type] - self.tel_types.append(type) + p_ped.append(ped_dict[tel_id][mask]) + p_spe.append(spe_dict[tel_id][mask]) + self.tel_id.append(tel_id) self.tel_pos_x[i] = tilt_coord[indices[i]].x.to(u.m).value @@ -705,7 +752,6 @@ class before minimisation can take place. This simply copies a ma.zeros(shape), ma.zeros(shape), ) - self.tel_types = np.array(self.tel_types) # Copy everything into our masked arrays for i in range(len(hillas_dict)): @@ -714,8 +760,8 @@ class before minimisation can take place. This simply copies a self.pixel_y[i][:array_len] = py[i] self.image[i][:array_len] = pa[i] self.time[i][:array_len] = pt[i] - self.ped[i][:] = self.ped_table[self.tel_types[i]] - self.spe[i][:] = 0.5 + self.ped[i][:array_len] = p_ped[i] + self.spe[i][:array_len] = p_spe[i] # Set the image mask mask = self.image <= 0.0 @@ -746,6 +792,8 @@ def predict( image_dict=None, mask_dict=None, time_dict=None, + ped_dict=None, + spe_dict=None, ): """Predict method for the ImPACT reconstructor. Used to calculate the reconstructed ImPACT shower geometry and energy. @@ -769,6 +817,8 @@ def predict( image_dict, time_dict, mask_dict, + ped_dict, + spe_dict, subarray, array_pointing, telescope_pointings, From 37c6bd4d2390a9df53553ff2bc1146c4cdb47bd9 Mon Sep 17 00:00:00 2001 From: gschwefer Date: Fri, 28 Feb 2025 14:00:45 +0100 Subject: [PATCH 04/43] Add tel_type key to dummy templates --- src/ctapipe/reco/impact_utilities.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/ctapipe/reco/impact_utilities.py b/src/ctapipe/reco/impact_utilities.py index e617c28807a..9e920521cc1 100644 --- a/src/ctapipe/reco/impact_utilities.py +++ b/src/ctapipe/reco/impact_utilities.py @@ -164,6 +164,7 @@ def generate_fake_template( def create_dummy_templates( output_file, energy, + tel_type, pe=1000.0, energy_range=np.logspace(-1, 1, 7), dist_range=np.linspace(0, 200, 5), @@ -192,5 +193,8 @@ def create_dummy_templates( template_dict[key] = template.T * pe + output_dict = {} + output_dict["data"] = template_dict + output_dict["tel_type"] = tel_type with gzip.open(output_file, "wb") as filehandler: - pickle.dump(template_dict, filehandler) + pickle.dump(output_dict, filehandler) From 5beb81df9b2fc3da896e143ea96d8c809f5ce866 Mon Sep 17 00:00:00 2001 From: gschwefer Date: Fri, 28 Feb 2025 14:01:16 +0100 Subject: [PATCH 05/43] Make tests pass gain --- src/ctapipe/reco/tests/test_ImPACT.py | 191 ++++++++++++++++++++++---- 1 file changed, 168 insertions(+), 23 deletions(-) diff --git a/src/ctapipe/reco/tests/test_ImPACT.py b/src/ctapipe/reco/tests/test_ImPACT.py index 3aa6f89d41b..cc0beb1817c 100644 --- a/src/ctapipe/reco/tests/test_ImPACT.py +++ b/src/ctapipe/reco/tests/test_ImPACT.py @@ -3,6 +3,7 @@ import pytest from astropy.coordinates import AltAz, Angle, SkyCoord from numpy.testing import assert_allclose +from traitlets.config import Config from ctapipe.containers import ( HillasParametersContainer, @@ -58,12 +59,45 @@ def setup_class(self): kurtosis=0, ) - def test_brightest_mean_average(self, example_subarray, table_profile): + def test_brightest_mean_average(self, tmp_path, example_subarray, table_profile): """ Test that averaging of the brightest pixel position give a sensible outcome """ - impact_reco = ImPACTReconstructor(example_subarray, table_profile) + for tel_type in example_subarray.telescope_types: + create_dummy_templates( + str(tmp_path) + "/{}.template.gz".format(str(tel_type)), + 1, + str(tel_type), + ) + + impact_config = Config( + { + "ImPACTReconstructor": { + "image_template_path": [ + [ + "type", + "MST_MST_FlashCam", + str(tmp_path) + "/MST_MST_FlashCam.template.gz", + ], + [ + "type", + "SST_ASTRI_ASTRICam", + str(tmp_path) + "/SST_ASTRI_ASTRICam.template.gz", + ], + [ + "type", + "LST_LST_LSTCam", + str(tmp_path) + "/LST_LST_LSTCam.template.gz", + ], + ] + } + } + ) + + impact_reco = ImPACTReconstructor( + example_subarray, table_profile, config=impact_config + ) pixel_x = np.array([0.0, 1.0, 0.0, -1.0]) * u.deg pixel_y = np.array([-1.0, 0.0, 1.0, 0.0]) * u.deg @@ -98,10 +132,43 @@ def test_translation(self): assert_allclose(xt, 1, rtol=0, atol=0.001) assert_allclose(yt, -1, rtol=0, atol=0.001) - def test_xmax_calculation(self, example_subarray, table_profile): + def test_xmax_calculation(self, tmp_path, example_subarray, table_profile): """Test calculation of hmax and interpolation of Xmax tables""" - impact_reco = ImPACTReconstructor(example_subarray, table_profile) + for tel_type in example_subarray.telescope_types: + create_dummy_templates( + str(tmp_path) + "/{}.template.gz".format(str(tel_type)), + 1, + str(tel_type), + ) + + impact_config = Config( + { + "ImPACTReconstructor": { + "image_template_path": [ + [ + "type", + "MST_MST_FlashCam", + str(tmp_path) + "/MST_MST_FlashCam.template.gz", + ], + [ + "type", + "SST_ASTRI_ASTRICam", + str(tmp_path) + "/SST_ASTRI_ASTRICam.template.gz", + ], + [ + "type", + "LST_LST_LSTCam", + str(tmp_path) + "/LST_LST_LSTCam.template.gz", + ], + ] + } + } + ) + + impact_reco = ImPACTReconstructor( + example_subarray, table_profile, config=impact_config + ) pixel_x = np.array([1, 1, 1]) * u.deg pixel_y = np.array([1, 1, 1]) * u.deg @@ -119,17 +186,45 @@ def test_xmax_calculation(self, example_subarray, table_profile): def test_interpolation(self, tmp_path, example_subarray, table_profile): """Test interpolation works on dummy template library""" - impact_reco = ImPACTReconstructor(example_subarray, table_profile) + for tel_type in example_subarray.telescope_types: + create_dummy_templates( + str(tmp_path) + "/{}.template.gz".format(str(tel_type)), + 1, + str(tel_type), + ) + + impact_config = Config( + { + "ImPACTReconstructor": { + "image_template_path": [ + [ + "type", + "MST_MST_FlashCam", + str(tmp_path) + "/MST_MST_FlashCam.template.gz", + ], + [ + "type", + "SST_ASTRI_ASTRICam", + str(tmp_path) + "/SST_ASTRI_ASTRICam.template.gz", + ], + [ + "type", + "LST_LST_LSTCam", + str(tmp_path) + "/LST_LST_LSTCam.template.gz", + ], + ] + } + } + ) + + impact_reco = ImPACTReconstructor( + example_subarray, table_profile, config=impact_config + ) - create_dummy_templates(str(tmp_path) + "/dummy.template.gz", 1) template, x, y = generate_fake_template(-1.5, 0.5) template *= 1000 - impact_reco.root_dir = str(tmp_path) - impact_reco.initialise_templates({1: "dummy"}) - - pred = impact_reco.image_prediction( - "dummy", + pred = impact_reco.prediction[(1, 2, 3, 4)]( 0, 0, np.array([1]), @@ -142,10 +237,40 @@ def test_interpolation(self, tmp_path, example_subarray, table_profile): assert_allclose(template.ravel() - pred, np.zeros_like(pred), atol=0.1) def test_fitting(self, tmp_path, example_subarray, table_profile): - impact_reco = ImPACTReconstructor(example_subarray, table_profile) + for tel_type in example_subarray.telescope_types: + create_dummy_templates( + str(tmp_path) + "/{}.template.gz".format(str(tel_type)), + 1, + str(tel_type), + ) + + impact_config = Config( + { + "ImPACTReconstructor": { + "image_template_path": [ + [ + "type", + "MST_MST_FlashCam", + str(tmp_path) + "/MST_MST_FlashCam.template.gz", + ], + [ + "type", + "SST_ASTRI_ASTRICam", + str(tmp_path) + "/SST_ASTRI_ASTRICam.template.gz", + ], + [ + "type", + "LST_LST_LSTCam", + str(tmp_path) + "/LST_LST_LSTCam.template.gz", + ], + ] + } + } + ) - create_dummy_templates(str(tmp_path) + "/dummy.template.gz", 1) - impact_reco.root_dir = str(tmp_path) + impact_reco = ImPACTReconstructor( + example_subarray, table_profile, config=impact_config + ) tel1, x, y = generate_fake_template(-1.5, 0.5, 0.3, 50, 50, ((-4, 4), (-4, 4))) tel2 = np.rot90(tel1) @@ -158,13 +283,16 @@ def test_fitting(self, tmp_path, example_subarray, table_profile): array_pointing = SkyCoord(alt=0 * u.deg, az=0 * u.deg, frame=AltAz) - impact_reco.tel_types = np.array(["dummy", "dummy", "dummy", "dummy"]) - impact_reco.initialise_templates( - {1: "dummy", 2: "dummy", 3: "dummy", 4: "dummy"} - ) + impact_reco.tel_id = [1, 2, 3, 4] + impact_reco.template_masks = { + tuple(np.arange(1, 5)): np.array([True, True, True, True]), + tuple(np.arange(5, 30)): np.array([False, False, False, False]), + tuple(np.arange(30, 99)): np.array([False, False, False, False]), + } impact_reco.zenith = 0 # *u.deg impact_reco.azimuth = 0 # *u.deg - impact_reco.ped = np.ones_like(image) # *u.deg + impact_reco.ped = np.ones_like(image) + impact_reco.spe = np.ones_like(image) impact_reco.image = image * 1000 impact_reco.hillas_parameters = [self.h1, self.h1, self.h1, self.h1] @@ -191,10 +319,27 @@ def test_selected_subarray( ): """test that reconstructor also works with "missing" ids""" - create_dummy_templates(str(tmp_path) + "/LSTCam.template.gz", 1) - subarray, event = subarray_and_event_gamma_off_axis_500_gev + create_dummy_templates( + str(tmp_path) + "/LST_LST_LSTCam.template.gz", + 1, + str(subarray.telescope_types[0]), + ) + impact_config = Config( + { + "ImPACTReconstructor": { + "image_template_path": [ + [ + "type", + "LST_LST_LSTCam", + str(tmp_path) + "/LST_LST_LSTCam.template.gz", + ] + ] + } + } + ) + shower_test = ReconstructedGeometryContainer() energy_test = ReconstructedEnergyContainer() shower_test.prefix = "test" @@ -214,8 +359,8 @@ def test_selected_subarray( event.dl2.stereo.geometry["test"] = shower_test event.dl2.stereo.energy["test_energy"] = energy_test - reconstructor = ImPACTReconstructor(subarray, table_profile) - reconstructor.root_dir = str(tmp_path) + reconstructor = ImPACTReconstructor(subarray, table_profile, config=impact_config) + reconstructor(event) assert event.dl2.stereo.geometry["ImPACTReconstructor"].is_valid assert event.dl2.stereo.energy["ImPACTReconstructor"].is_valid From 8fbb6fcf5b1326ad85e4914d624be6f16a67e765 Mon Sep 17 00:00:00 2001 From: gschwefer Date: Fri, 28 Feb 2025 14:02:54 +0100 Subject: [PATCH 06/43] Remove dummy reconstructor option --- src/ctapipe/reco/impact.py | 11 +---------- 1 file changed, 1 insertion(+), 10 deletions(-) diff --git a/src/ctapipe/reco/impact.py b/src/ctapipe/reco/impact.py index 13200d699bb..cf717760508 100644 --- a/src/ctapipe/reco/impact.py +++ b/src/ctapipe/reco/impact.py @@ -77,7 +77,6 @@ "SST_SST_SST-Camera": 0.5, "SST_SST_CHEC": 0.5, "SST_ASTRI_ASTRICam": 0.5, - "dummy": 0.01, "UNKNOWN-960PX": 1.0, } @@ -88,7 +87,6 @@ "SST_SST_SST-Camera": 0.6, "SST_SST_CHEC": 0.6, "SST_ASTRI_ASTRICam": 0.6, - "dummy": 0.6, "UNKNOWN-960PX": 0.6, } @@ -120,9 +118,6 @@ class ImPACTReconstructor(HillasGeometryReconstructor): The telescope subarray to use for reconstruction atmosphere_profile : ctapipe.atmosphere.AtmosphereDensityProfile Density vs. altitude profile of the local atmosphere - dummy_reconstructor : bool, optional - Option to use a set of dummy templates. This can be used for testing the algorithm, - but for any actual reconstruction should be set to its default False References ---------- @@ -162,9 +157,7 @@ class ImPACTReconstructor(HillasGeometryReconstructor): property = ReconstructionProperty.ENERGY | ReconstructionProperty.GEOMETRY - def __init__( - self, subarray, atmosphere_profile, dummy_reconstructor=False, **kwargs - ): + def __init__(self, subarray, atmosphere_profile, **kwargs): if Minuit is None: raise OptionalDependencyMissing("iminuit") from None @@ -204,8 +197,6 @@ def __init__( self.array_direction = None self.nominal_frame = None - self.dummy_reconstructor = dummy_reconstructor - def __call__(self, event): """ Perform the full shower geometry reconstruction on the input event. From c4bdb2b0a9c9a354dc8c9102df459d326e324829 Mon Sep 17 00:00:00 2001 From: gschwefer Date: Mon, 3 Mar 2025 15:51:33 +0100 Subject: [PATCH 07/43] Update shower processor config to work with ImPACT --- .../reco/tests/test_shower_processor.py | 37 ++++++++++++++++++- 1 file changed, 35 insertions(+), 2 deletions(-) diff --git a/src/ctapipe/reco/tests/test_shower_processor.py b/src/ctapipe/reco/tests/test_shower_processor.py index 1560bf837e4..58b14f3098c 100644 --- a/src/ctapipe/reco/tests/test_shower_processor.py +++ b/src/ctapipe/reco/tests/test_shower_processor.py @@ -10,6 +10,7 @@ from ctapipe.calib import CameraCalibrator from ctapipe.image import ImageProcessor from ctapipe.reco import ShowerProcessor +from ctapipe.reco.impact_utilities import create_dummy_templates from ctapipe.utils import get_dataset_path SIMTEL_PATH = get_dataset_path( @@ -42,13 +43,44 @@ def table_profile(): ids=["HillasIntersection", "HillasReconstructor", "ImPACTReconstructor"], ) def test_shower_processor_geometry( - example_event, example_subarray, reconstructor_types, table_profile + example_event, example_subarray, reconstructor_types, table_profile, tmp_path ): """Ensure we get shower geometry when we input an event with parametrized images.""" calibrate = CameraCalibrator(subarray=example_subarray) - config = Config() + for tel_type in example_subarray.telescope_types: + create_dummy_templates( + str(tmp_path) + "/{}.template.gz".format(str(tel_type)), + 1, + str(tel_type), + ) + + config = Config( + { + "ShowerProcessor": { + "ImPACTReconstructor": { + "image_template_path": [ + [ + "type", + "MST_MST_FlashCam", + str(tmp_path) + "/MST_MST_FlashCam.template.gz", + ], + [ + "type", + "SST_ASTRI_ASTRICam", + str(tmp_path) + "/SST_ASTRI_ASTRICam.template.gz", + ], + [ + "type", + "LST_LST_LSTCam", + str(tmp_path) + "/LST_LST_LSTCam.template.gz", + ], + ] + } + } + } + ) process_images = ImageProcessor( subarray=example_subarray, image_cleaner_type="MARSImageCleaner" @@ -58,6 +90,7 @@ def test_shower_processor_geometry( subarray=example_subarray, atmosphere_profile=table_profile, reconstructor_types=reconstructor_types, + config=config, ) calibrate(example_event) From f87bf317bce5b7e160cba158aa785018fc84463a Mon Sep 17 00:00:00 2001 From: gschwefer Date: Mon, 3 Mar 2025 16:08:26 +0100 Subject: [PATCH 08/43] Add template paths to reconstruction tests --- .../reco/tests/test_reconstruction_methods.py | 40 ++++++++++++++++--- 1 file changed, 34 insertions(+), 6 deletions(-) diff --git a/src/ctapipe/reco/tests/test_reconstruction_methods.py b/src/ctapipe/reco/tests/test_reconstruction_methods.py index ba80f17b208..1f8c8cefef1 100644 --- a/src/ctapipe/reco/tests/test_reconstruction_methods.py +++ b/src/ctapipe/reco/tests/test_reconstruction_methods.py @@ -1,18 +1,20 @@ import pytest from astropy import units as u +from traitlets.config import Config from ctapipe.calib import CameraCalibrator from ctapipe.image import ImageProcessor from ctapipe.io import EventSource from ctapipe.reco import HillasIntersection, HillasReconstructor from ctapipe.reco.impact import ImPACTReconstructor +from ctapipe.reco.impact_utilities import create_dummy_templates from ctapipe.utils import get_dataset_path reconstructors = [HillasIntersection, HillasReconstructor, ImPACTReconstructor] @pytest.mark.parametrize("cls", reconstructors) -def test_reconstructors(cls): +def test_reconstructors(cls, tmp_path): """ a test of the complete fit procedure on one event including: • tailcut cleaning @@ -24,7 +26,6 @@ def test_reconstructors(cls): filename = get_dataset_path( "gamma_LaPalma_baseline_20Zd_180Az_prod3b_test.simtel.gz" ) - template_file = get_dataset_path("LSTCam.template.gz") source = EventSource(filename, max_events=4, focal_length_choice="EQUIVALENT") subarray = source.subarray @@ -35,15 +36,42 @@ def test_reconstructors(cls): calib(event) image_processor(event) + reco_config = Config() + if cls is ImPACTReconstructor: pytest.importorskip("iminuit") + for tel_type in subarray.telescope_types: + create_dummy_templates( + str(tmp_path) + "/{}.template.gz".format(str(tel_type)), + 1, + str(tel_type), + ) + + reco_config = Config( + { + "ImPACTReconstructor": { + "image_template_path": [ + [ + "type", + "LST_LST_LSTCam", + str(tmp_path) + "/LST_LST_LSTCam.template.gz", + ], + [ + "type", + "MST_MST_NectarCam", + str(tmp_path) + "/MST_MST_NectarCam.template.gz", + ], + ] + } + } + ) + reconstructor = cls( - subarray, atmosphere_profile=source.atmosphere_density_profile + subarray, + atmosphere_profile=source.atmosphere_density_profile, + config=reco_config, ) - if cls is ImPACTReconstructor: - reconstructor.root_dir = str(template_file.parents[0]) - reconstructor(event) name = cls.__name__ From 17bb07f8233ccec7db4f9f55bf150c41d62ba80e Mon Sep 17 00:00:00 2001 From: gschwefer Date: Mon, 3 Mar 2025 16:23:08 +0100 Subject: [PATCH 09/43] Add changelog --- docs/changes/2705.feature.rst | 1 + 1 file changed, 1 insertion(+) create mode 100644 docs/changes/2705.feature.rst diff --git a/docs/changes/2705.feature.rst b/docs/changes/2705.feature.rst new file mode 100644 index 00000000000..1bf70b68862 --- /dev/null +++ b/docs/changes/2705.feature.rst @@ -0,0 +1 @@ +Make template paths and likelihood parameters TelecopeParameters in ImPACT From f37c75b1611e4d9e4c0a252048975cd2224715d7 Mon Sep 17 00:00:00 2001 From: gschwefer Date: Fri, 21 Mar 2025 09:40:26 +0100 Subject: [PATCH 10/43] Make template bounds a global variable --- src/ctapipe/reco/impact.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/ctapipe/reco/impact.py b/src/ctapipe/reco/impact.py index cf717760508..0ce0b900a53 100644 --- a/src/ctapipe/reco/impact.py +++ b/src/ctapipe/reco/impact.py @@ -90,6 +90,10 @@ "UNKNOWN-960PX": 0.6, } +# Validity bounds of the templates in degree in camera coordinates. +# Needed to initialize the templates +TEMPLATE_BOUNDS = ((-5, 1), (-1.5, 1.5)) + __all__ = ["ImPACTReconstructor"] @@ -332,7 +336,7 @@ def set_up_templates(self): for template_path, tel_ids in template_sort_dict.items(): net_interpolator = TemplateNetworkInterpolator( - template_path, bounds=((-5, 1), (-1.5, 1.5)) + template_path, bounds=TEMPLATE_BOUNDS ) interp_tel_string = net_interpolator.tel_type_string From a2542cc0c52ea3af7d69d98adb58a5ead7a8c4fb Mon Sep 17 00:00:00 2001 From: Dan Parsons Date: Thu, 6 Nov 2025 14:20:42 +0100 Subject: [PATCH 11/43] Commit minimal version of FreePACT code --- src/ctapipe/reco/__init__.py | 3 + src/ctapipe/reco/freepact.py | 158 ++++++ src/ctapipe/reco/impact.py | 175 ++----- src/ctapipe/reco/impact_utilities.py | 26 +- .../utils/template_network_interpolator.py | 463 +++++++++++++++--- 5 files changed, 617 insertions(+), 208 deletions(-) create mode 100644 src/ctapipe/reco/freepact.py diff --git a/src/ctapipe/reco/__init__.py b/src/ctapipe/reco/__init__.py index 267b68b4b5c..ca90b69f447 100644 --- a/src/ctapipe/reco/__init__.py +++ b/src/ctapipe/reco/__init__.py @@ -2,6 +2,7 @@ # reconstructors must be imported before ShowerProcessor, so # they are available there +from .freepact import FreePACTProtonReconstructor, FreePACTReconstructor from .hillas_intersection import HillasIntersection from .hillas_reconstructor import HillasReconstructor from .impact import ImPACTReconstructor @@ -26,6 +27,8 @@ "ShowerProcessor", "HillasReconstructor", "ImPACTReconstructor", + "FreePACTReconstructor", + "FreePACTProtonReconstructor", "HillasIntersection", "EnergyRegressor", "ParticleClassifier", diff --git a/src/ctapipe/reco/freepact.py b/src/ctapipe/reco/freepact.py new file mode 100644 index 00000000000..2819a44492f --- /dev/null +++ b/src/ctapipe/reco/freepact.py @@ -0,0 +1,158 @@ +"""FreePACT Reconstructor for ctapipe +This module implements the FreePACT reconstructor, which is a subclass of ImPACTReconstructor. +It uses a neural network to predict th likelihood camera images based on the shower parameters.""" + +import numpy as np +import numpy.ma as ma + +from ctapipe.core import traits +from ctapipe.core.telescope_component import TelescopeParameter +from ctapipe.reco.impact import ImPACTReconstructor +from ctapipe.utils.template_network_interpolator import FreePACTInterpolator + +from ..compat import COPY_IF_NEEDED +from .impact_utilities import rotate_translate + +__all__ = ["FreePACTReconstructor"] + + +class FreePACTReconstructor(ImPACTReconstructor): + """ + FreePACT Reconstructor + This class implements the FreePACT reconstructor, which uses a neural network + to predict the likelihood of camera images based on shower parameters. + """ + + image_template_path = TelescopeParameter( + trait=traits.Path(exists=True, directory_ok=True, allow_none=False), + allow_none=False, + help=("Path to the image templates to be used in the reconstruction"), + ).tag(config=True) + + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + + def set_up_templates(self): + """Set up the templates for the FreePACT reconstructor.""" + template_sort_dict = {} + + for tel_id in self.subarray.tel_ids: + if self.image_template_path.tel[tel_id] not in template_sort_dict.keys(): + template_sort_dict[self.image_template_path.tel[tel_id]] = [tel_id] + else: + template_sort_dict[self.image_template_path.tel[tel_id]].append(tel_id) + + for template_path, tel_ids in template_sort_dict.items(): + net_interpolator = FreePACTInterpolator(template_path) + + self.prediction[tuple(tel_ids)] = net_interpolator + + def get_likelihood( + self, + source_x, + source_y, + core_x, + core_y, + energy, + x_max_scale, + goodness_of_fit=False, + ): + """Get the likelihood that the image predicted at the given test + position matches the camera image. + + Parameters + ---------- + source_x: float + Source position of shower in the nominal system (in deg) + source_y: float + Source position of shower in the nominal system (in deg) + core_x: float + Core position of shower in tilted telescope system (in m) + core_y: float + Core position of shower in tilted telescope system (in m) + energy: float + Shower energy (in TeV) + x_max_scale: float + Scaling factor applied to geometrically calculated Xmax + + Returns + ------- + float: Likelihood the model represents the camera image at this position + + """ + if np.isnan(source_x) or np.isnan(source_y): + return 1e8 + + # First we add units back onto everything. Currently not + # handled very well, maybe in future we could just put + # everything in the correct units when loading in the class + # and ignore them from then on + + zenith = self.zenith + azimuth = self.azimuth + + # Geometrically calculate the depth of maximum given this test position + x_max_guess = self.get_shower_max(source_x, source_y, core_x, core_y, zenith) + x_max_guess *= x_max_scale + + # Convert to binning of Xmax + x_max_diff = x_max_guess + + # Calculate impact distance for all telescopes + impact = np.sqrt( + (self.tel_pos_x - core_x) ** 2 + (self.tel_pos_y - core_y) ** 2 + ) + # And the expected rotation angle + phi = np.arctan2(self.tel_pos_y - core_y, self.tel_pos_x - core_x) + + # Rotate and translate all pixels such that they match the + # template orientation + # numba does not support masked arrays, work on underlying array and add mask back + pix_x_rot, pix_y_rot = rotate_translate( + self.pixel_y.data, self.pixel_x.data, source_y, source_x, -phi + ) + pix_x_rot = np.ma.array(pix_x_rot, mask=self.pixel_x.mask, copy=COPY_IF_NEEDED) + pix_y_rot = np.ma.array(pix_y_rot, mask=self.pixel_y.mask, copy=COPY_IF_NEEDED) + + # In the interpolator class we can gain speed advantages by using masked arrays + # so we need to make sure here everything is masked + likelihood = ma.zeros(self.image.shape) + likelihood.mask = ma.getmask(self.image) + + for tel_ids, template in self.prediction.items(): + template_mask = self.template_masks[tel_ids] + if np.any(template_mask): + likelihood[template_mask] = template( + np.rad2deg(zenith), + azimuth, + energy * np.ones_like(impact[template_mask]), + impact[template_mask], + x_max_diff * np.ones_like(impact[template_mask]), + np.rad2deg(pix_x_rot[template_mask]), + np.rad2deg(pix_y_rot[template_mask]), + self.image[template_mask], + ) + if goodness_of_fit: + # return -2 * np.sum(likelihood[self.image>5]) / np.sum(self.image>5) + # print(ma.getmask(self.image)) + return -2 * ma.sum(likelihood) / np.sum(~ma.getmask(self.image)) + + likelihood = ma.sum(likelihood) * -2 + + return likelihood + + +class FreePACTProtonReconstructor(FreePACTReconstructor): + """FreePACT Proton Reconstructor + This class implements the FreePACT reconstructor for proton showers. + It uses a neural network to predict the likelihood of camera images based on shower parameters. + """ + + image_template_path = TelescopeParameter( + trait=traits.Path(exists=True, directory_ok=True, allow_none=False), + allow_none=False, + help=("Path to the image templates to be used in the reconstruction"), + ).tag(config=True) + + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) diff --git a/src/ctapipe/reco/impact.py b/src/ctapipe/reco/impact.py index 0ce0b900a53..955833016ef 100644 --- a/src/ctapipe/reco/impact.py +++ b/src/ctapipe/reco/impact.py @@ -8,7 +8,6 @@ import numpy.ma as ma from astropy import units as u from astropy.coordinates import AltAz, SkyCoord -from scipy.stats import norm from ctapipe.core import traits from ctapipe.core.telescope_component import TelescopeParameter @@ -24,16 +23,12 @@ project_to_ground, ) from ..core import Provenance -from ..fitting import lts_linear_regression from ..image.cleaning import dilate from ..image.pixel_likelihood import ( mean_poisson_likelihood_gaussian, neg_log_likelihood_approx, ) -from ..utils.template_network_interpolator import ( - TemplateNetworkInterpolator, - TimeGradientInterpolator, -) +from ..utils.template_network_interpolator import TemplateNetworkInterpolator from .impact_utilities import ( EmptyImages, create_seed, @@ -54,16 +49,6 @@ PROV = Provenance() -INVALID_GEOMETRY = ReconstructedGeometryContainer( - telescopes=[], - prefix="ImPACTReconstructor", -) - -INVALID_ENERGY = ReconstructedEnergyContainer( - prefix="ImPACTReconstructor", - telescopes=[], -) - # These are settings for the iminuit minimizer MINUIT_ERRORDEF = 0.5 # 0.5 for a log-likelihood cost function for correct errors MINUIT_STRATEGY = 1 # Default minimization strategy, 2 is careful, 0 is fast @@ -77,6 +62,9 @@ "SST_SST_SST-Camera": 0.5, "SST_SST_CHEC": 0.5, "SST_ASTRI_ASTRICam": 0.5, + "MST_HESS-1_HESS-1U": 1.0, + "LST_H.E.S.S. CT5 (876 mirrors)_H.E.S.S. CT5 with FlashCam": 1.0, + "dummy": 0.01, "UNKNOWN-960PX": 1.0, } @@ -87,7 +75,10 @@ "SST_SST_SST-Camera": 0.6, "SST_SST_CHEC": 0.6, "SST_ASTRI_ASTRICam": 0.6, + "MST_HESS-1_HESS-1U": 0.6, + "LST_H.E.S.S. CT5 (876 mirrors)_H.E.S.S. CT5 with FlashCam": 0.6, "UNKNOWN-960PX": 0.6, + "dummy": 0.6, } # Validity bounds of the templates in degree in camera coordinates. @@ -135,12 +126,14 @@ class ImPACTReconstructor(HillasGeometryReconstructor): help=("Path to the image templates to be used in the reconstruction"), ).tag(config=True) - # The time gradient templates are optional, so None is allowed here - time_gradient_template_path = TelescopeParameter( - trait=traits.Path(exists=True, directory_ok=False, allow_none=True), + # The SPE and pedestal width parameters are also configurable as TelescopeParameters. + # None is allowed and the default value. In that case, either values from the event monitoring data are used + # or if that is also not available, a value from a hardcoded backup dict is used. + + pedestal_width = traits.FloatTelescopeParameter( allow_none=True, default_value=None, - help=("Path to the time gradient templates to be used in the reconstruction"), + help="Pedestal width parameter for the likelihood", ).tag(config=True) # The SPE and pedestal width parameters are also configurable as TelescopeParameters. @@ -200,6 +193,17 @@ def __init__(self, subarray, atmosphere_profile, **kwargs): self.array_direction = None self.nominal_frame = None + self.site_height = self.subarray.reference_location.height.to_value(u.m) + + self.INVALID_GEOMETRY = ReconstructedGeometryContainer( + telescopes=[], + prefix=self.__class__.__name__, + ) + + self.INVALID_ENERGY = ReconstructedEnergyContainer( + prefix=self.__class__.__name__, + telescopes=[], + ) def __call__(self, event): """ @@ -214,8 +218,8 @@ def __call__(self, event): try: hillas_dict = self._create_hillas_dict(event) except (TooFewTelescopesException, InvalidWidthException): - event.dl2.stereo.geometry[self.__class__.__name__] = INVALID_GEOMETRY - event.dl2.stereo.energy[self.__class__.__name__] = INVALID_ENERGY + event.dl2.stereo.geometry[self.__class__.__name__] = self.INVALID_GEOMETRY + event.dl2.stereo.energy[self.__class__.__name__] = self.INVALID_ENERGY self._store_impact_parameter(event) return @@ -258,16 +262,14 @@ def __call__(self, event): ) # Dilate the images around the original cleaning to help the fit - for _ in range(3): + for _ in range(4): mask = dilate(self.subarray.tel[tel_id].camera.geometry, mask) mask_dict[tel_id] = mask # Next, we look for geometry and energy seeds from previously applied reconstructors. # Both need to be present at elast once for ImPACT to run. - reco_geom_pred = event.dl2.stereo.geometry - valid_geometry_seed = False for geom_pred in reco_geom_pred.values(): if geom_pred.is_valid: @@ -275,7 +277,6 @@ def __call__(self, event): break reco_energy_pred = event.dl2.stereo.energy - valid_energy_seed = False for E_pred in reco_energy_pred.values(): if E_pred.is_valid: @@ -283,8 +284,8 @@ def __call__(self, event): break if valid_geometry_seed is False or valid_energy_seed is False: - event.dl2.stereo.geometry[self.__class__.__name__] = INVALID_GEOMETRY - event.dl2.stereo.energy[self.__class__.__name__] = INVALID_ENERGY + event.dl2.stereo.geometry[self.__class__.__name__] = self.INVALID_GEOMETRY + event.dl2.stereo.energy[self.__class__.__name__] = self.INVALID_ENERGY self._store_impact_parameter(event) return @@ -307,10 +308,8 @@ def __call__(self, event): self._store_impact_parameter(event) def set_up_templates(self): + """Set up the templates for the ImPACT reconstruction.""" template_sort_dict = {} - time_template_sort_dict = {} - - self.use_time_gradient = True for tel_id in self.subarray.tel_ids: if self.image_template_path.tel[tel_id] not in template_sort_dict.keys(): @@ -318,22 +317,6 @@ def set_up_templates(self): else: template_sort_dict[self.image_template_path.tel[tel_id]].append(tel_id) - if self.time_gradient_template_path.tel[tel_id] is not None: - if ( - self.time_gradient_template_path.tel[tel_id] - not in time_template_sort_dict.keys() - ): - time_template_sort_dict[ - self.time_gradient_template_path.tel[tel_id] - ] = [tel_id] - else: - time_template_sort_dict[ - self.time_gradient_template_path.tel[tel_id] - ].append(tel_id) - - else: - self.use_time_gradient = False - for template_path, tel_ids in template_sort_dict.items(): net_interpolator = TemplateNetworkInterpolator( template_path, bounds=TEMPLATE_BOUNDS @@ -348,20 +331,6 @@ def set_up_templates(self): self.prediction[tuple(tel_ids)] = net_interpolator - if self.use_time_gradient: - for template_path, tel_ids in time_template_sort_dict.items(): - time_interpolator = TimeGradientInterpolator(template_path) - - interp_tel_string = time_interpolator.tel_type_string - - for id in tel_ids: - if interp_tel_string != str(self.subarray.tel[id]): - raise ValueError( - "You are using templates that are not intended for this telescope type" - ) - - self.time_prediction[tuple(tel_ids)] = time_interpolator - def get_hillas_mean(self): """This is a simple function to find the peak position of each image in an event which will be used later in the Xmax calculation. Peak is @@ -436,10 +405,7 @@ def get_shower_max(self, source_x, source_y, core_x, core_y, zen): ) / np.sum(weight) # Add on the height of the detector above sea level - mean_height_asl = ( - mean_height_above_ground - + self.subarray.reference_location.height.to_value(u.m) - ) + mean_height_asl = mean_height_above_ground + self.site_height if mean_height_asl > 100000 or np.isnan(mean_height_asl): mean_height_asl = 100000 @@ -448,7 +414,6 @@ def get_shower_max(self, source_x, source_y, core_x, core_y, zen): slant_depth = self.atmosphere_profile.slant_depth_from_height( mean_height_asl * u.m, zen * u.rad ) - return slant_depth.to_value(u.g / (u.cm * u.cm)) def get_likelihood( @@ -533,10 +498,6 @@ def get_likelihood( prediction = ma.zeros(self.image.shape) prediction.mask = ma.getmask(self.image) - time_gradients, time_gradients_uncertainty = ( - np.zeros(self.image.shape[0]), - np.zeros(self.image.shape[0]), - ) # Loop over all telescope types and get prediction for tel_ids, template in self.prediction.items(): @@ -552,49 +513,6 @@ def get_likelihood( np.rad2deg(pix_y_rot[template_mask]), ) - if self.use_time_gradient: - for tel_ids, time_template in self.time_prediction.items(): - time_template_mask = self.time_template_masks[tel_ids] - if np.any(time_template_mask): - time_pred = time_template( - np.rad2deg(zenith), - azimuth, - energy * np.ones_like(impact[time_template_mask]), - impact[time_template_mask], - x_max_diff * np.ones_like(impact[time_template_mask]), - ) - - time_gradients[time_template_mask] = time_pred.T[0] - time_gradients_uncertainty[time_template_mask] = time_pred.T[1] - - time_gradients_uncertainty[time_gradients_uncertainty == 0] = 1e-6 - - chi2 = 0 - for telescope_index, (image, time) in enumerate(zip(self.image, self.time)): - time_mask = np.logical_and( - np.invert(ma.getmask(image)), - time > 0, - ) - time_mask = np.logical_and(time_mask, np.isfinite(time)) - time_mask = np.logical_and(time_mask, image > 5) - if ( - np.sum(time_mask) > 3 - and time_gradients_uncertainty[telescope_index] > 0 - ): - time_slope = lts_linear_regression( - x=np.rad2deg(pix_x_rot[telescope_index][time_mask]), - y=time[time_mask], - samples=3, - )[0][0] - - time_like = -1 * norm.logpdf( - time_slope, - loc=time_gradients[telescope_index], - scale=time_gradients_uncertainty[telescope_index], - ) - - chi2 += time_like - # Likelihood function will break if we find a NaN or a 0 prediction[np.isnan(prediction)] = 1e-8 prediction[prediction < 1e-8] = 1e-8 @@ -620,8 +538,6 @@ def get_likelihood( like = np.sum(like) final_sum = like - if self.use_time_gradient: - final_sum += chi2 return final_sum @@ -768,11 +684,7 @@ class before minimisation can take place. This simply copies a tels_key: np.isin(list(hillas_dict.keys()), tels_key) for tels_key in self.prediction.keys() } - if self.use_time_gradient: - self.time_template_masks = { - tels_key: np.isin(list(hillas_dict.keys()), tels_key) - for tels_key in self.time_prediction.keys() - } + # Finally run some functions to get ready for the event self.get_hillas_mean() @@ -827,6 +739,8 @@ def predict( for cleaning in shower_seed: # Copy all of our seed parameters out of the shower objects # We need to convert the shower direction to the nominal system + if cleaning == "FreePACTReconstructor": + continue horizon_seed = SkyCoord( az=shower_seed[cleaning].az, alt=shower_seed[cleaning].alt, @@ -851,8 +765,14 @@ def predict( zenith = 90 * u.deg - self.array_direction.alt for energy_reco in energy_seed: + # We need to convert the shower direction to the nominal system + + if energy_reco == "FreePACTReconstructor": + continue energy = energy_seed[energy_reco].energy.value + # for i in range(1): + # energy = 0.5 # TEST!!! seed, step, limits = create_seed( source_x, source_y, tilt_x, tilt_y, energy ) @@ -869,18 +789,18 @@ def predict( like_min = like if fit_params is None: - return INVALID_GEOMETRY, INVALID_ENERGY + return self.INVALID_GEOMETRY, self.INVALID_ENERGY # Now do full minimisation seed = create_seed( fit_params[0], fit_params[1], fit_params[2], fit_params[3], fit_params[4] ) - fit_params, errors, like = self.minimise( - params=seed[0], - step=seed[1], - limits=seed[2], - ) + # fit_params, errors, like = self.minimise( + # params=seed[0], + # step=seed[1], + # limits=seed[2], + # ) # Create a container class for reconstructed shower @@ -986,6 +906,7 @@ def minimise( ------- tuple: best fit parameters and errors """ + limits = np.asarray(limits) energy = params[4] @@ -998,7 +919,7 @@ def minimise( source_y=params[1], core_x=params[2], core_y=params[3], - energy=energy, + energy=energy / 0.8, x_max_scale=xmax_scale, goodness_of_fit=False, ) diff --git a/src/ctapipe/reco/impact_utilities.py b/src/ctapipe/reco/impact_utilities.py index 9e920521cc1..3577d7282a6 100644 --- a/src/ctapipe/reco/impact_utilities.py +++ b/src/ctapipe/reco/impact_utilities.py @@ -22,7 +22,7 @@ class EmptyImages(Exception): def guess_shower_depth(energy): """ Simple estimation of depth of shower max based on the expected gamma-ray elongation - rate. + rate Parameters ---------- @@ -67,16 +67,16 @@ def rotate_translate(pixel_pos_x, pixel_pos_y, x_trans, y_trans, phi): pixel_pos_trans_x, pixel_pos_trans_y = np.zeros(shape), np.zeros(shape) for i in range(shape[0]): - cosine_angle = np.cos(phi[i]) - sin_angle = np.sin(phi[i]) - for j in range(shape[1]): - pixel_pos_trans_x[i][j] = (x_trans - pixel_pos_x[i][j]) * cosine_angle - ( - y_trans - pixel_pos_y[i][j] - ) * sin_angle - pixel_pos_trans_y[i][j] = (pixel_pos_x[i][j] - x_trans) * sin_angle + ( - pixel_pos_y[i][j] - y_trans - ) * cosine_angle + cosine_angle = np.cos(phi[i][j]) + sin_angle = np.sin(phi[i][j]) + + pixel_pos_trans_x[i][j] = ( + x_trans[i][j] - pixel_pos_x[i][j] + ) * cosine_angle - (y_trans[i][j] - pixel_pos_y[i][j]) * sin_angle + pixel_pos_trans_y[i][j] = ( + pixel_pos_x[i][j] - x_trans[i][j] + ) * sin_angle + (pixel_pos_y[i][j] - y_trans[i][j]) * cosine_angle return pixel_pos_trans_x, pixel_pos_trans_y @@ -116,15 +116,15 @@ def create_seed(source_x, source_y, tilt_x, tilt_y, energy): seed = [source_x, source_y, tilt_x, tilt_y, en_seed, 1.2] # Take a reasonable first guess at step size - step = [0.04 / 57.3, 0.04 / 57.3, 10, 10, en_seed * 0.05, 0.05, 0.01] + step = [0.05 / 57.3, 0.05 / 57.3, 10, 10, en_seed * 0.1, 0.1, 0.01] # And some sensible limits of the fit range limits = [ [source_x - 1.5 / 57.3, source_x + 1.5 / 57.3], [source_y - 1.5 / 57.3, source_y + 1.5 / 57.3], [tilt_x - 100, tilt_x + 100], [tilt_y - 100, tilt_y + 100], - [lower_en_limit, en_seed * 2], - [0.8, 1.2], + [lower_en_limit, en_seed * 10], + [0.5, 1.5], [0.0, 0.01], ] diff --git a/src/ctapipe/utils/template_network_interpolator.py b/src/ctapipe/utils/template_network_interpolator.py index 78578af06c5..e1477d3af9c 100644 --- a/src/ctapipe/utils/template_network_interpolator.py +++ b/src/ctapipe/utils/template_network_interpolator.py @@ -1,11 +1,17 @@ import gzip +import os import pickle +import re +import numba import numpy as np import numpy.ma as ma +import tensorflow as tf from .unstructured_interpolator import UnstructuredInterpolator +tf.config.set_visible_devices([], "GPU") + class BaseTemplate: """ @@ -93,28 +99,24 @@ def _get_bounds(self, zenith, azimuth): # If we are not at the edge we need to reduce the index by one if index != 0: - index -= 1 index_upper -= 1 - zenith_bound = (index, index_upper) - # Do the same again for azimuth angle if len(self.azimuths) == 0: azimuth_bound = self.azimuths[0], self.azimuths[0] else: - index = np.searchsorted(self.azimuths, azimuth) + az_index = np.searchsorted(self.azimuths, azimuth) # Except in this case we need to loop back around rather than stay at the edge - if index != 0: - index -= 1 + if az_index != 0: + az_index -= 1 - index_upper = index - if index != len(self.azimuths) - 1: - index_upper += 1 + az_index_upper = az_index + if az_index != len(self.azimuths) - 1: + az_index_upper += 1 else: - index_upper = 0 - - azimuth_bound = (index, index_upper) + az_index_upper = 0 + azimuth_bound = (az_index, az_index_upper) # Return our boundaries return zenith_bound, azimuth_bound @@ -157,74 +159,49 @@ def _create_interpolator(self, zenith_bin, azimuth_bin): def perform_interpolation(self, zenith, azimuth, interpolation_array, points): """ - Parameters - ---------- - zenith - azimuth - interpolation_array - points - Returns - ------- + Interpolate values for given zenith and azimuth using bilinear interpolation. """ - zenith_bounds, azimuth_bounds = self._get_bounds(zenith, azimuth) + zl, zu = zenith_bounds + al, au = azimuth_bounds - zenith_lower, zenith_upper = zenith_bounds - azimuth_lower, azimuth_upper = azimuth_bounds + def eval(z, a): + # Evaluate or create the interpolator for the given zenith and azimuth bin + return self._evaluate_interpolator(z, a, interpolation_array, points) - # First lower azimuth bound - if self.interpolator[zenith_lower][azimuth_lower] is None: - self._create_interpolator(zenith_lower, azimuth_lower) - evaluate_azimuth_lower1 = self.interpolator[zenith_lower][azimuth_lower]( - interpolation_array, points - ) + v_ll = eval(zl, al) + if len(self.zeniths) == 1 and len(self.azimuths) == 1: + # If only one zenith and azimuth, return directly + return v_ll - if len(self.zeniths) == 1 and len(self.zeniths) == 1: - return evaluate_azimuth_lower1 - - if self.interpolator[zenith_upper][azimuth_lower] is None: - self._create_interpolator(zenith_upper, azimuth_lower) - evaluate_azimuth_lower2 = self.interpolator[zenith_upper][azimuth_lower]( - interpolation_array, points - ) + v_ul = eval(zu, al) + v_lu = eval(zl, au) + v_uu = eval(zu, au) - evaluate_azimuth_lower = self._linear_interpolation( - zenith, - self.zeniths[zenith_lower], - self.zeniths[zenith_upper], - evaluate_azimuth_lower1, - evaluate_azimuth_lower2, + # Interpolate along zenith for lower and upper azimuth bounds + v_l = self._linear_interpolation( + zenith, self.zeniths[zl], self.zeniths[zu], v_ll, v_ul ) - # Then the upper - if self.interpolator[zenith_lower][azimuth_upper] is None: - self._create_interpolator(zenith_lower, azimuth_upper) - evaluate_azimuth_upper1 = self.interpolator[zenith_lower][azimuth_upper]( - interpolation_array, points + v_u = self._linear_interpolation( + zenith, self.zeniths[zl], self.zeniths[zu], v_lu, v_uu ) - if self.interpolator[zenith_upper][azimuth_upper] is None: - self._create_interpolator(zenith_upper, azimuth_upper) - evaluate_azimuth_upper2 = self.interpolator[zenith_upper][azimuth_upper]( - interpolation_array, points + # Interpolate along azimuth between the two zenith-interpolated values + return self._linear_interpolation( + azimuth, self.azimuths[al], self.azimuths[au], v_l, v_u ) - evaluate_azimuth_upper = self._linear_interpolation( - zenith, - self.zeniths[zenith_lower], - self.zeniths[zenith_upper], - evaluate_azimuth_upper1, - evaluate_azimuth_upper2, - ) - # And finally interpolate between the azimuths - result = self._linear_interpolation( - azimuth, - self.azimuths[azimuth_lower], - self.azimuths[azimuth_upper], - evaluate_azimuth_lower, - evaluate_azimuth_upper, - ) + def _evaluate_interpolator( + self, zenith_bin, azimuth_bin, interpolation_array, points + ): + """ + Helper function to create and evaluate the interpolator for a given zenith and azimuth bin. + If the interpolator does not exist, it is created. + """ + if self.interpolator[zenith_bin][azimuth_bin] is None: + self._create_interpolator(zenith_bin, azimuth_bin) - return result + return self.interpolator[zenith_bin][azimuth_bin](interpolation_array, points) @staticmethod def _linear_interpolation(point, grid1, grid2, value1, value2): @@ -414,17 +391,367 @@ def __eq__(self, other): return self.template_file == other.template_file +def load_prediction_files_filtered(directory): + """ + Finds all files in the given directory matching the pattern: + predict_{telescope}_{zenith}deg_{azimuth}azm_{offset}off_{species}.keras + Only includes files where telescope and species match the specified values. + + Returns a dictionary where the key is a tuple of (zenith, azimuth, offset) + and the value is the absolute file path. + """ + pattern = re.compile( + r"predict_(?P[^_]+)_(?P\d+)deg_(?P\d+)deg_(?P[\d\.]+)off.keras" + ) + result = {} + for filename in os.listdir(directory): + match = pattern.match(filename) + if match: + key = ( + int(match.group("zenith")), + int(match.group("azimuth")), + float(match.group("offset")), + ) + abs_path = os.path.abspath(os.path.join(directory, filename)) + model = tf.keras.models.load_model(abs_path) + model.layers[-1].activation = tf.keras.activations.linear + result[key] = model + return result + + +@numba.njit +def custom_symlog(value, linear_threshold=10.0): + """ + Apply a symmetric logarithm transformation to the input array. This is + implemented using numba due to potential performance constraints. + + Parameters + ---------- + value: ndarray + Input array to transform + linear_threshold: float + Threshold below which to apply linear scaling + """ + for i in range(value.shape[0]): + for j in range(value.shape[1]): + if np.abs(value[i][j]) < linear_threshold: + value[i][j] = value[i][j] + else: + value[i][j] = np.sign(value[i][j]) * ( + np.log2(np.abs(value[i][j] / linear_threshold)) + linear_threshold + ) + + return value + + +@tf.function +def evaluate_model(model, parameters): + """ + Evaluate the model with the given parameters. TF decorator used to + execute so model is compiled into a callable TensorFlow graph. + + Parameters + ---------- + model: tf.keras.Model + The model to evaluate + parameters: ndarray + The parameters to pass to the model + + Returns + ------- + tf.Tensor + The output of the model + """ + return model(parameters) + + +@tf.function +def evaluate_model_interpolate( + model_ul, + model_uu, + model_lu, + model_ll, + zenith, + azimuth, + zenith_u, + zenith_l, + azimuth_u, + azimuth_l, + parameters, +): + """ + Evaluate the model with the given parameters in an interpolated grid. TF decorator used to + execute so model is compiled into a callable TensorFlow graph. + + Parameters + ---------- + model_ul: tf.keras.Model + The model to evaluate for the upper left corner + model_uu: tf.keras.Model + The model to evaluate for the upper right corner + model_lu: tf.keras.Model + The model to evaluate for the lower right corner + model_ll: tf.keras.Model + The model to evaluate for the lower left corner + zenith: float + The zenith angle to evaluate + azimuth: float + The azimuth angle to evaluate + zenith_u: float + The upper zenith bound for interpolation + zenith_l: float + The lower zenith bound for interpolation + azimuth_u: float + The upper azimuth bound for interpolation + azimuth_l: float + The lower azimuth bound for interpolation + parameters: ndarray + The parameters to pass to the model + + Returns + ------- + tf.Tensor + The output of the model + """ + v_ll = model_ll(parameters) + v_ul = model_ul(parameters) + v_lu = model_lu(parameters) + v_uu = model_uu(parameters) + # print(v_ll, zenith) + v_l = ((v_ul - v_ll) * (zenith - zenith_l) / (zenith_u - zenith_l)) + v_ll + v_u = ((v_uu - v_lu) * (zenith - zenith_l) / (zenith_u - zenith_l)) + v_lu + + # Interpolate along azimuth between the two zenith-interpolated values + value = ((v_u - v_l) * (azimuth - azimuth_l) / (azimuth_u - azimuth_l)) + v_l + + return value + + +class FreePACTInterpolator(BaseTemplate): + """ + Class for interpolating between the FreePACT predictions + """ + + def __init__(self, directory): + """ + Parameters + ---------- + directory: str + Directory containing the FreePACT prediction files + """ + + super().__init__() + + data_input_dict = load_prediction_files_filtered(directory) + # self.tel_type_string = telescope_type + + keys = np.array(list(data_input_dict.keys())) + values = list(data_input_dict.values()) + + self.no_zenaz = False + + # First check if we even have a zen and azimuth entry + if len(keys) > 1: + # If we do then for the sake of speed lets + self._create_table_matrix(keys, values) + else: + # If not we work as before + # Currently impact is not set up for offset dependent templates. + # Therefore remove offset (last) dimension from interpolator + self.interpolator = values[0] + self.no_zenaz = True + + def _create_interpolator(self, zenith_bin, azimuth_bin): + """ + Creates unstructured interpolator object in a given zenith and azimuth bin when + needed + Parameters + ---------- + zenith_bin: int + Bin number of requested zenith + azimuth_bin: int + Bin number of requested azimuth + Returns + ------- + None + """ + # Get our requested zenith and azimuth + + zenith = self.zeniths[zenith_bin] + azimuth = self.azimuths[azimuth_bin] + + # Select these values from our range of keys + selection = np.logical_and(self.keys.T[0] == zenith, self.keys.T[1] == azimuth) + bin_sel = np.where(selection)[0][0] + + # Create interpolator using this selection + # Currently freepact is not set up for offset dependent templates. + # Therefore remove offset (last) dimension from interpolator + self.interpolator[zenith_bin][azimuth_bin] = self.values[bin_sel] + + # We can now remove these entries. + self.keys = self.keys[np.invert(selection)] + del self.values[bin_sel] + # self.values = self.values[np.invert(selection)] + + def perform_interpolation(self, zenith, azimuth, interpolation_array, points): + """ + Perform interpolation between template models for a given zenith and azimuth angle. + + Parameters + ---------- + zenith: float + Zenith angle in degrees + azimuth: float + Azimuth angle in degrees + interpolation_array: array-like + Empty in this case + points: array-like + Array of points to evaluate + + Returns + ------- + array-like + Interpolated values + """ + zenith_bounds, azimuth_bounds = self._get_bounds(zenith, azimuth) + zl, zu = zenith_bounds + al, au = azimuth_bounds + + if self.interpolator[zl][al] is None: + self._create_interpolator(zl, al) + if self.interpolator[zu][au] is None: + self._create_interpolator(zu, au) + if self.interpolator[zu][al] is None: + self._create_interpolator(zu, al) + if self.interpolator[zl][au] is None: + self._create_interpolator(zl, au) + + value = evaluate_model_interpolate( + self.interpolator[zu][al], + self.interpolator[zu][au], + self.interpolator[zl][au], + self.interpolator[zl][al], + tf.convert_to_tensor( + zenith * np.ones((points.shape[0], 1)), dtype=tf.float32 + ), + tf.convert_to_tensor( + azimuth * np.ones((points.shape[0], 1)), dtype=tf.float32 + ), + tf.convert_to_tensor( + self.zeniths[zu] * np.ones((points.shape[0], 1)), + dtype=tf.float32, + ), + tf.convert_to_tensor( + self.zeniths[zl] * np.ones((points.shape[0], 1)), + dtype=tf.float32, + ), + tf.convert_to_tensor( + self.azimuths[au] * np.ones((points.shape[0], 1)), + dtype=tf.float32, + ), + tf.convert_to_tensor( + self.azimuths[al] * np.ones((points.shape[0], 1)), + dtype=tf.float32, + ), + tf.convert_to_tensor(points), + ) + + return value.numpy() # np.asarray(value) + + def _evaluate_interpolator( + self, zenith_bin, azimuth_bin, interpolation_array, points + ): + """ + Helper function to create and evaluate the interpolator for a given zenith and azimuth bin. + If the interpolator does not exist, it is created. + + Parameters + ---------- + zenith_bin: int + Zenith bin index + azimuth_bin: int + Azimuth bin index + interpolation_array: array-like + Array of points to interpolate + points: array-like + Array of points to evaluate + """ + if self.interpolator[zenith_bin][azimuth_bin] is None: + self._create_interpolator(zenith_bin, azimuth_bin) + return np.asarray( + evaluate_model(self.interpolator[zenith_bin][azimuth_bin], points) + ) + + # return self.interpolator[zenith_bin][azimuth_bin].predict(interpolation_array, verbose=0, batch_size=10000) + + def __call__(self, zenith, azimuth, energy, impact, xmax, xb, yb, amplitude): + """ + Evaluate expected time gradient for a set of shower parameters and pixel positions + Parameters + ---------- + energy: array-like + Energy of interpolated template + impact: array-like + Impact distance of interpolated template + xmax: array-like + Depth of maximum of interpolated templates + Returns + ------- + ndarray: Time Gradient expectation and RMS values + """ + shape = xb.shape + repeat_num = xb.shape[-1] + + amplitude = custom_symlog(ma.getdata(amplitude) / 10, linear_threshold=2) + impact = impact / 100 + energy = np.log10(energy) + xmax = xmax / 100 + + array = np.stack( + ( + xb.ravel(), + yb.ravel(), + np.repeat(impact, repeat_num, axis=-1), + np.repeat(energy, repeat_num, axis=-1), + np.repeat(xmax, repeat_num, axis=-1), + amplitude.ravel(), + ), + axis=-1, + ) + + if self.no_zenaz: + interpolated_value = evaluate_model(self.interpolator, array).numpy() + interpolated_value = interpolated_value.reshape(shape) + else: + interpolated_value = self.perform_interpolation( + zenith, azimuth, None, array + ) + interpolated_value = interpolated_value.reshape(shape) + + return interpolated_value + + def reset(self): + return True + + class DummyTemplateInterpolator: + """Dummy template interpolator for testing purposes.""" + def __call__(self, zenith, azimuth, energy, impact, xmax, xb, yb): return np.ones_like(xb) def reset(self): + """Reset the interpolator.""" return True class DummyTimeInterpolator: + """Dummy time interpolator for testing purposes.""" + def __call__(self, energy, impact, xmax): return np.ones_like(energy) def reset(self): + """Reset the interpolator.""" return True From c03acdcac5da0a9f109df3561aa240daf4bb9d1e Mon Sep 17 00:00:00 2001 From: Dan Parsons Date: Fri, 7 Nov 2025 11:52:19 +0100 Subject: [PATCH 12/43] Minor commit to change to new container structure. --- src/ctapipe/reco/impact.py | 15 +++++++++++++-- src/ctapipe/reco/impact_utilities.py | 17 ++++++++--------- 2 files changed, 21 insertions(+), 11 deletions(-) diff --git a/src/ctapipe/reco/impact.py b/src/ctapipe/reco/impact.py index 955833016ef..a18c0236312 100644 --- a/src/ctapipe/reco/impact.py +++ b/src/ctapipe/reco/impact.py @@ -241,12 +241,23 @@ def __call__(self, event): time_dict[tel_id] = event.dl1.tel[tel_id].peak_time mask = event.dl1.tel[tel_id].image_mask if self.pedestal_width.tel[tel_id] is None: - if event.mon.tel[tel_id].pedestal.charge_std is None: + selected_gain_channel = event.dl0.tel[tel_id].selected_gain_channel + + if ( + event.monitoring.tel[ + tel_id + ].camera.pixel_statistics.sky_pedestal_image.std + is None + ): ped_dict[tel_id] = BACKUP_PED_TABLE[ str(self.subarray.tel[tel_id]) ] * np.ones(self.subarray.tel[tel_id].camera.geometry.n_pixels) else: - ped_dict[tel_id] = event.mon.tel[tel_id].pedestal.charge_std[0] + ped_dict[tel_id] = event.monitoring.tel[ + tel_id + ].camera.pixel_statistics.sky_pedestal_image.std[ + selected_gain_channel, : + ] else: ped_dict[tel_id] = self.pedestal_width.tel[tel_id] * np.ones( self.subarray.tel[tel_id].camera.geometry.n_pixels diff --git a/src/ctapipe/reco/impact_utilities.py b/src/ctapipe/reco/impact_utilities.py index 3577d7282a6..7e96736564b 100644 --- a/src/ctapipe/reco/impact_utilities.py +++ b/src/ctapipe/reco/impact_utilities.py @@ -67,16 +67,15 @@ def rotate_translate(pixel_pos_x, pixel_pos_y, x_trans, y_trans, phi): pixel_pos_trans_x, pixel_pos_trans_y = np.zeros(shape), np.zeros(shape) for i in range(shape[0]): + cosine_angle = np.cos(phi[i]) + sin_angle = np.sin(phi[i]) for j in range(shape[1]): - cosine_angle = np.cos(phi[i][j]) - sin_angle = np.sin(phi[i][j]) - - pixel_pos_trans_x[i][j] = ( - x_trans[i][j] - pixel_pos_x[i][j] - ) * cosine_angle - (y_trans[i][j] - pixel_pos_y[i][j]) * sin_angle - pixel_pos_trans_y[i][j] = ( - pixel_pos_x[i][j] - x_trans[i][j] - ) * sin_angle + (pixel_pos_y[i][j] - y_trans[i][j]) * cosine_angle + pixel_pos_trans_x[i][j] = (x_trans - pixel_pos_x[i][j]) * cosine_angle - ( + y_trans - pixel_pos_y[i][j] + ) * sin_angle + pixel_pos_trans_y[i][j] = (pixel_pos_x[i][j] - x_trans) * sin_angle + ( + pixel_pos_y[i][j] - y_trans + ) * cosine_angle return pixel_pos_trans_x, pixel_pos_trans_y From 4514e786127ff93a69c764d0823015da6678dfb8 Mon Sep 17 00:00:00 2001 From: Dan Parsons Date: Wed, 21 Jan 2026 09:54:34 +0100 Subject: [PATCH 13/43] Add tensorflow as optional dependency --- pyproject.toml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index d50f546f9e8..b204ecb9b1b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -58,7 +58,8 @@ all = [ "ctapipe[eventio]", "iminuit >=2", "matplotlib ~=3.0", - "pyirf ~=0.13.0" + "pyirf ~=0.13.0", + "tensorflow >=2.19.0" ] tests = [ From 4b34c72742091fca3d7a782b10117033b573930a Mon Sep 17 00:00:00 2001 From: Dan Parsons Date: Wed, 21 Jan 2026 11:07:34 +0100 Subject: [PATCH 14/43] Fix dependency versions in pyproject.toml --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index b204ecb9b1b..61ca8b4086a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -59,7 +59,7 @@ all = [ "iminuit >=2", "matplotlib ~=3.0", "pyirf ~=0.13.0", - "tensorflow >=2.19.0" + "tensorflow >=2.15" ] tests = [ From c2226a33c5006833d042b5cd567f7691bfeb09a3 Mon Sep 17 00:00:00 2001 From: Dan Parsons Date: Thu, 22 Jan 2026 13:20:09 +0100 Subject: [PATCH 15/43] Fixing CI build --- .github/workflows/ci.yml | 2 +- pyproject.toml | 39 ++++++++++++++------------------------- 2 files changed, 15 insertions(+), 26 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 0396381a611..b2415d0a376 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -58,7 +58,7 @@ jobs: python-version: "3.12" install-method: pip extra-args: ["codecov"] - extras: tests,all + extras: tests,all, impact - name: Linux (3.13, mamba) os: ubuntu-24.04 diff --git a/pyproject.toml b/pyproject.toml index 61ca8b4086a..17c0506fee9 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -17,29 +17,30 @@ license = "BSD-3-Clause" classifiers = [ "Intended Audience :: Science/Research", "Programming Language :: Python :: 3", + "Programming Language :: Python :: 3.10", + "Programming Language :: Python :: 3.11", "Programming Language :: Python :: 3.12", "Programming Language :: Python :: 3.13", - "Programming Language :: Python :: 3.14", "Programming Language :: Python :: Implementation :: CPython", "Topic :: Scientific/Engineering :: Astronomy", "Development Status :: 3 - Alpha", ] dynamic = ["version"] -requires-python = ">=3.12" +requires-python = ">=3.10" dependencies = [ - "astropy >=6.1,<8.0.0", + "astropy >=6.1,<8.0.0a0", "docutils", "joblib", - "numba >=0.59.0", - "numpy >=1.26,<3.0.0", + "numba >=0.57", + "numpy >=1.24,<3.0.0a0", "packaging", "psutil", "pyyaml >=5.1", "requests", "scikit-learn !=1.4.0", # 1.4.0 breaks with astropy tables, before and after works - "scipy ~=1.15", + "scipy ~=1.10", "tables ~=3.4", "tqdm >=4.32", "traitlets ~=5.6", @@ -47,25 +48,20 @@ dependencies = [ [project.optional-dependencies] -eventio = [ - "eventio >=1.9.1,<3.0.0", -] - # all is with all optional *runtime* dependencies # use `dev` to get really all dependencies all = [ "bokeh ~=3.0", - "ctapipe[eventio]", + "eventio >=1.9.1,<2.0.0a0", "iminuit >=2", "matplotlib ~=3.0", "pyirf ~=0.13.0", - "tensorflow >=2.15" ] tests = [ # at the moment, essentially all tests rely on test data from simtel # it doesn't make sense to skip all of these. - "ctapipe[eventio]", + "eventio >=1.9.1,<2.0.0a0", "h5py", "pandas", "pytest >= 7.0", @@ -76,6 +72,9 @@ tests = [ "astroquery", ] +impact = ["tensorflow >=2.19.0", + "iminuit >=2"] + docs = [ "ctapipe[all]", "ffmpeg-python", @@ -94,6 +93,7 @@ docs = [ "sphinx_automodapi", "sphinxcontrib-bibtex", "sphinx-changelog", + "tomli; python_version < '3.11'", ] dev = [ @@ -162,19 +162,8 @@ filterwarnings = [ "error::ctapipe.utils.deprecation.CTAPipeDeprecationWarning", "error::ctapipe.instrument.FromNameWarning", "error::ctapipe.core.traits.NoneDefaultNotAllowedWarning", - "error::tables.UnclosedFileWarning", "ignore:`np.MachAr` is deprecated:DeprecationWarning", "ignore::ctapipe.core.provenance.MissingReferenceMetadata", - # inside matplotlib, already fixed but occurs in "oldest-deps" tests - "ignore:'parseString':DeprecationWarning:matplotlib", - "ignore:'resetCache':DeprecationWarning:matplotlib", - "ignore:'oneOf':DeprecationWarning:matplotlib", - "ignore:'leaveWhiteSpace':DeprecationWarning:matplotlib", - "ignore:'setName':DeprecationWarning:matplotlib", - "ignore:'setParseAction':DeprecationWarning:matplotlib", - "ignore:'endQuoteChar':DeprecationWarning:matplotlib", - "ignore:'unquoteResults':DeprecationWarning:matplotlib", - "ignore:'parseAll':DeprecationWarning:pyparsing", ] norecursedirs = [ ".git", @@ -226,7 +215,7 @@ omit = ["src/ctapipe/_version.py"] [tool.ruff] -target-version = "py312" +target-version = "py310" line-length = 88 [tool.ruff.lint] From 0250307ead8e9d6a9799002ea542dc23f2a7fcba Mon Sep 17 00:00:00 2001 From: Dan Parsons Date: Thu, 22 Jan 2026 13:20:59 +0100 Subject: [PATCH 16/43] Reformat of file --- .github/workflows/ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index b2415d0a376..0cdc3148e8a 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -58,7 +58,7 @@ jobs: python-version: "3.12" install-method: pip extra-args: ["codecov"] - extras: tests,all, impact + extras: tests,all,impact - name: Linux (3.13, mamba) os: ubuntu-24.04 From d039bf3fe518cd32ab0c56a9769a71f7c6a54498 Mon Sep 17 00:00:00 2001 From: Dan Parsons Date: Thu, 22 Jan 2026 13:27:02 +0100 Subject: [PATCH 17/43] Fixing CI build... again --- .github/workflows/ci.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 0cdc3148e8a..7c6591bbc00 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -44,14 +44,14 @@ jobs: os: ubuntu-24.04 python-version: "3.12" install-method: mamba - extras: tests,all + extras: tests,all,impact extra-args: ["oldest-deps"] - name: Linux (3.12, pip, minimal-dependencies) os: ubuntu-24.04 python-version: "3.12" install-method: pip - extras: tests + extras: tests,impact - name: Linux (3.12, pip, coverage) os: ubuntu-24.04 From ee0e95b5f8a73737e425e5872f28573817ee946f Mon Sep 17 00:00:00 2001 From: Dan Parsons Date: Fri, 23 Jan 2026 09:44:51 +0100 Subject: [PATCH 18/43] Making import of tensorflow optional --- .github/workflows/ci.yml | 2 +- pyproject.toml | 1 - src/ctapipe/utils/template_network_interpolator.py | 11 ++++++++++- 3 files changed, 11 insertions(+), 3 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 7c6591bbc00..1c0e555f886 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -76,7 +76,7 @@ jobs: os: macos-15-intel python-version: "3.12" install-method: pip - extras: tests,all + extras: tests,all,impact - name: macos (3.14, arm64, mamba) os: macos-15 diff --git a/pyproject.toml b/pyproject.toml index 17c0506fee9..3b945d98b09 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -53,7 +53,6 @@ dependencies = [ all = [ "bokeh ~=3.0", "eventio >=1.9.1,<2.0.0a0", - "iminuit >=2", "matplotlib ~=3.0", "pyirf ~=0.13.0", ] diff --git a/src/ctapipe/utils/template_network_interpolator.py b/src/ctapipe/utils/template_network_interpolator.py index e1477d3af9c..ebb104b6bf5 100644 --- a/src/ctapipe/utils/template_network_interpolator.py +++ b/src/ctapipe/utils/template_network_interpolator.py @@ -6,7 +6,13 @@ import numba import numpy as np import numpy.ma as ma -import tensorflow as tf + +from ..exceptions import OptionalDependencyMissing + +try: + import tensorflow as tf +except ModuleNotFoundError: + tf = None from .unstructured_interpolator import UnstructuredInterpolator @@ -561,6 +567,9 @@ def __init__(self, directory): self.interpolator = values[0] self.no_zenaz = True + if tf is None: + raise OptionalDependencyMissing("tensorflow") + def _create_interpolator(self, zenith_bin, azimuth_bin): """ Creates unstructured interpolator object in a given zenith and azimuth bin when From 64b92f2e10e8bbe79bed7a2649f06ce233b63826 Mon Sep 17 00:00:00 2001 From: Dan Parsons Date: Fri, 23 Jan 2026 09:57:52 +0100 Subject: [PATCH 19/43] Fix messed up merge --- pyproject.toml | 42 +++++++++++++++++++++++++++--------------- 1 file changed, 27 insertions(+), 15 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 3b945d98b09..5ed3b594b5b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -17,30 +17,29 @@ license = "BSD-3-Clause" classifiers = [ "Intended Audience :: Science/Research", "Programming Language :: Python :: 3", - "Programming Language :: Python :: 3.10", - "Programming Language :: Python :: 3.11", "Programming Language :: Python :: 3.12", "Programming Language :: Python :: 3.13", + "Programming Language :: Python :: 3.14", "Programming Language :: Python :: Implementation :: CPython", "Topic :: Scientific/Engineering :: Astronomy", "Development Status :: 3 - Alpha", ] dynamic = ["version"] -requires-python = ">=3.10" +requires-python = ">=3.12" dependencies = [ - "astropy >=6.1,<8.0.0a0", + "astropy >=6.1,<8.0.0", "docutils", "joblib", - "numba >=0.57", - "numpy >=1.24,<3.0.0a0", + "numba >=0.59.0", + "numpy >=1.26,<3.0.0", "packaging", "psutil", "pyyaml >=5.1", "requests", "scikit-learn !=1.4.0", # 1.4.0 breaks with astropy tables, before and after works - "scipy ~=1.10", + "scipy ~=1.15", "tables ~=3.4", "tqdm >=4.32", "traitlets ~=5.6", @@ -48,19 +47,23 @@ dependencies = [ [project.optional-dependencies] +eventio = [ + "eventio >=1.9.1,<3.0.0", +] + # all is with all optional *runtime* dependencies # use `dev` to get really all dependencies all = [ "bokeh ~=3.0", - "eventio >=1.9.1,<2.0.0a0", + "ctapipe[eventio]", "matplotlib ~=3.0", - "pyirf ~=0.13.0", + "pyirf ~=0.13.0" ] tests = [ # at the moment, essentially all tests rely on test data from simtel # it doesn't make sense to skip all of these. - "eventio >=1.9.1,<2.0.0a0", + "ctapipe[eventio]", "h5py", "pandas", "pytest >= 7.0", @@ -71,9 +74,6 @@ tests = [ "astroquery", ] -impact = ["tensorflow >=2.19.0", - "iminuit >=2"] - docs = [ "ctapipe[all]", "ffmpeg-python", @@ -92,7 +92,6 @@ docs = [ "sphinx_automodapi", "sphinxcontrib-bibtex", "sphinx-changelog", - "tomli; python_version < '3.11'", ] dev = [ @@ -101,6 +100,8 @@ dev = [ "setuptools_scm[toml]", ] +impact = ["tensorflow >=2.16.0", + "iminuit >=2"] [project.scripts] ctapipe-info = "ctapipe.tools.info:main" @@ -161,8 +162,19 @@ filterwarnings = [ "error::ctapipe.utils.deprecation.CTAPipeDeprecationWarning", "error::ctapipe.instrument.FromNameWarning", "error::ctapipe.core.traits.NoneDefaultNotAllowedWarning", + "error::tables.UnclosedFileWarning", "ignore:`np.MachAr` is deprecated:DeprecationWarning", "ignore::ctapipe.core.provenance.MissingReferenceMetadata", + # inside matplotlib, already fixed but occurs in "oldest-deps" tests + "ignore:'parseString':DeprecationWarning:matplotlib", + "ignore:'resetCache':DeprecationWarning:matplotlib", + "ignore:'oneOf':DeprecationWarning:matplotlib", + "ignore:'leaveWhiteSpace':DeprecationWarning:matplotlib", + "ignore:'setName':DeprecationWarning:matplotlib", + "ignore:'setParseAction':DeprecationWarning:matplotlib", + "ignore:'endQuoteChar':DeprecationWarning:matplotlib", + "ignore:'unquoteResults':DeprecationWarning:matplotlib", + "ignore:'parseAll':DeprecationWarning:pyparsing", ] norecursedirs = [ ".git", @@ -214,7 +226,7 @@ omit = ["src/ctapipe/_version.py"] [tool.ruff] -target-version = "py310" +target-version = "py312" line-length = 88 [tool.ruff.lint] From ffd8d550a2e8c1c551e7fb36684fe82d213e1ec6 Mon Sep 17 00:00:00 2001 From: Dan Parsons Date: Fri, 23 Jan 2026 10:09:50 +0100 Subject: [PATCH 20/43] Fixing optional dependency --- src/ctapipe/utils/template_network_interpolator.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/ctapipe/utils/template_network_interpolator.py b/src/ctapipe/utils/template_network_interpolator.py index ebb104b6bf5..4e2ae4c9a25 100644 --- a/src/ctapipe/utils/template_network_interpolator.py +++ b/src/ctapipe/utils/template_network_interpolator.py @@ -12,7 +12,7 @@ try: import tensorflow as tf except ModuleNotFoundError: - tf = None + raise OptionalDependencyMissing("tensorflow") from .unstructured_interpolator import UnstructuredInterpolator @@ -547,6 +547,8 @@ def __init__(self, directory): """ super().__init__() + if tf is None: + raise OptionalDependencyMissing("tensorflow") data_input_dict = load_prediction_files_filtered(directory) # self.tel_type_string = telescope_type @@ -567,9 +569,6 @@ def __init__(self, directory): self.interpolator = values[0] self.no_zenaz = True - if tf is None: - raise OptionalDependencyMissing("tensorflow") - def _create_interpolator(self, zenith_bin, azimuth_bin): """ Creates unstructured interpolator object in a given zenith and azimuth bin when From 35b22adf1ae6308357e25cb5259c99a2dd77fbce Mon Sep 17 00:00:00 2001 From: Dan Parsons Date: Fri, 23 Jan 2026 11:24:09 +0100 Subject: [PATCH 21/43] Fixed failing unit tests --- src/ctapipe/reco/freepact.py | 8 ++++++-- src/ctapipe/reco/impact.py | 6 ++++-- src/ctapipe/utils/template_network_interpolator.py | 8 ++++---- 3 files changed, 14 insertions(+), 8 deletions(-) diff --git a/src/ctapipe/reco/freepact.py b/src/ctapipe/reco/freepact.py index 2819a44492f..d40f7344064 100644 --- a/src/ctapipe/reco/freepact.py +++ b/src/ctapipe/reco/freepact.py @@ -24,7 +24,9 @@ class FreePACTReconstructor(ImPACTReconstructor): """ image_template_path = TelescopeParameter( - trait=traits.Path(exists=True, directory_ok=True, allow_none=False), + trait=traits.Path( + exists=True, directory_ok=True, allow_none=False, default_value="./" + ), allow_none=False, help=("Path to the image templates to be used in the reconstruction"), ).tag(config=True) @@ -149,7 +151,9 @@ class FreePACTProtonReconstructor(FreePACTReconstructor): """ image_template_path = TelescopeParameter( - trait=traits.Path(exists=True, directory_ok=True, allow_none=False), + trait=traits.Path( + exists=True, directory_ok=True, allow_none=False, default_value="./" + ), allow_none=False, help=("Path to the image templates to be used in the reconstruction"), ).tag(config=True) diff --git a/src/ctapipe/reco/impact.py b/src/ctapipe/reco/impact.py index 6f7095946f9..4cfe5dce576 100644 --- a/src/ctapipe/reco/impact.py +++ b/src/ctapipe/reco/impact.py @@ -122,7 +122,9 @@ class ImPACTReconstructor(HillasGeometryReconstructor): """ image_template_path = TelescopeParameter( - trait=traits.Path(exists=True, directory_ok=False, allow_none=False), + trait=traits.Path( + exists=True, directory_ok=False, allow_none=False, default_value="./" + ), allow_none=False, help=("Path to the image templates to be used in the reconstruction"), ).tag(config=True) @@ -247,7 +249,7 @@ def __call__(self, event): if ( event.monitoring.tel[ tel_id - ].camera.pixel_statistics.sky_pedestal_image.std + ].camera.pixel_statistics.pedestal_image.std is None ): ped_dict[tel_id] = BACKUP_PED_TABLE[ diff --git a/src/ctapipe/utils/template_network_interpolator.py b/src/ctapipe/utils/template_network_interpolator.py index 4e2ae4c9a25..43eaedd0946 100644 --- a/src/ctapipe/utils/template_network_interpolator.py +++ b/src/ctapipe/utils/template_network_interpolator.py @@ -92,8 +92,8 @@ def _get_bounds(self, zenith, azimuth): # If we only have one zenith angle available in the templates, don't bother # searching - if len(self.zeniths) == 0: - zenith_bound = self.zeniths[0], self.zeniths[0] + if len(self.zeniths) == 1: + zenith_bound = (0, 0) else: # Otherwise search for where our zenith lies in the available range index = np.searchsorted(self.zeniths, zenith, side="left") @@ -108,8 +108,8 @@ def _get_bounds(self, zenith, azimuth): index_upper -= 1 zenith_bound = (index, index_upper) # Do the same again for azimuth angle - if len(self.azimuths) == 0: - azimuth_bound = self.azimuths[0], self.azimuths[0] + if len(self.azimuths) == 1: + azimuth_bound = (0, 0) else: az_index = np.searchsorted(self.azimuths, azimuth) # Except in this case we need to loop back around rather than stay at the edge From 2f468ecf82b65bbf14903863f092aa552f9ebe20 Mon Sep 17 00:00:00 2001 From: Dan Parsons Date: Fri, 23 Jan 2026 14:56:49 +0100 Subject: [PATCH 22/43] Removed pedestal and spe default table --- src/ctapipe/reco/impact.py | 79 ++++++--------------------- src/ctapipe/reco/tests/test_ImPACT.py | 4 +- 2 files changed, 19 insertions(+), 64 deletions(-) diff --git a/src/ctapipe/reco/impact.py b/src/ctapipe/reco/impact.py index 4cfe5dce576..0c6641b1132 100644 --- a/src/ctapipe/reco/impact.py +++ b/src/ctapipe/reco/impact.py @@ -56,32 +56,6 @@ MINUIT_TOLERANCE_FACTOR = 1000 # Tolerance for convergence according to EDM criterion MIGRAD_ITERATE = 1 # Do not call migrad again if convergence was not reached -BACKUP_PED_TABLE = { - "LST_LST_LSTCam": 1.4, - "MST_MST_NectarCam": 1.3, - "MST_MST_FlashCam": 1.3, - "SST_SST_SST-Camera": 0.5, - "SST_SST_CHEC": 0.5, - "SST_ASTRI_ASTRICam": 0.5, - "MST_HESS-1_HESS-1U": 1.0, - "LST_H.E.S.S. CT5 (876 mirrors)_H.E.S.S. CT5 with FlashCam": 1.0, - "dummy": 0.01, - "UNKNOWN-960PX": 1.0, -} - -BACKUP_SPE_TABLE = { - "LST_LST_LSTCam": 0.6, - "MST_MST_NectarCam": 0.6, - "MST_MST_FlashCam": 0.6, - "SST_SST_SST-Camera": 0.6, - "SST_SST_CHEC": 0.6, - "SST_ASTRI_ASTRICam": 0.6, - "MST_HESS-1_HESS-1U": 0.6, - "LST_H.E.S.S. CT5 (876 mirrors)_H.E.S.S. CT5 with FlashCam": 0.6, - "UNKNOWN-960PX": 0.6, - "dummy": 0.6, -} - # Validity bounds of the templates in degree in camera coordinates. # Needed to initialize the templates TEMPLATE_BOUNDS = ((-5, 1), (-1.5, 1.5)) @@ -139,16 +113,6 @@ class ImPACTReconstructor(HillasGeometryReconstructor): help="Pedestal width parameter for the likelihood", ).tag(config=True) - # The SPE and pedestal width parameters are also configurable as TelescopeParameters. - # None is allowed and the default value. In that case, either values from the event monitoring data are used - # or if that is also not available, a value from a hardcoded backup dict is used. - - pedestal_width = traits.FloatTelescopeParameter( - allow_none=True, - default_value=None, - help="Pedestal width parameter for the likelihood", - ).tag(config=True) - spe_width = traits.FloatTelescopeParameter( allow_none=True, default_value=None, @@ -243,37 +207,26 @@ def __call__(self, event): image_dict[tel_id] = image time_dict[tel_id] = event.dl1.tel[tel_id].peak_time mask = event.dl1.tel[tel_id].image_mask - if self.pedestal_width.tel[tel_id] is None: - selected_gain_channel = event.dl0.tel[tel_id].selected_gain_channel - - if ( - event.monitoring.tel[ - tel_id - ].camera.pixel_statistics.pedestal_image.std - is None - ): - ped_dict[tel_id] = BACKUP_PED_TABLE[ - str(self.subarray.tel[tel_id]) - ] * np.ones(self.subarray.tel[tel_id].camera.geometry.n_pixels) - else: - ped_dict[tel_id] = event.monitoring.tel[ - tel_id - ].camera.pixel_statistics.sky_pedestal_image.std[ - selected_gain_channel, : - ] - else: + + selected_gain_channel = event.dl0.tel[tel_id].selected_gain_channel + + if ( + event.monitoring.tel[tel_id].camera.pixel_statistics.pedestal_image.std + is None + ): ped_dict[tel_id] = self.pedestal_width.tel[tel_id] * np.ones( self.subarray.tel[tel_id].camera.geometry.n_pixels ) - - if self.spe_width.tel[tel_id] is None: - spe_dict[tel_id] = BACKUP_SPE_TABLE[ - str(self.subarray.tel[tel_id]) - ] * np.ones(self.subarray.tel[tel_id].camera.geometry.n_pixels) else: - spe_dict[tel_id] = self.spe_width.tel[tel_id] * np.ones( - self.subarray.tel[tel_id].camera.geometry.n_pixels - ) + ped_dict[tel_id] = event.monitoring.tel[ + tel_id + ].camera.pixel_statistics.sky_pedestal_image.std[ + selected_gain_channel, : + ] + + spe_dict[tel_id] = self.spe_width.tel[tel_id] * np.ones( + self.subarray.tel[tel_id].camera.geometry.n_pixels + ) # Dilate the images around the original cleaning to help the fit for _ in range(4): diff --git a/src/ctapipe/reco/tests/test_ImPACT.py b/src/ctapipe/reco/tests/test_ImPACT.py index cc0beb1817c..448b4ea63dc 100644 --- a/src/ctapipe/reco/tests/test_ImPACT.py +++ b/src/ctapipe/reco/tests/test_ImPACT.py @@ -335,7 +335,9 @@ def test_selected_subarray( "LST_LST_LSTCam", str(tmp_path) + "/LST_LST_LSTCam.template.gz", ] - ] + ], + "pedestal_width": [["type", "LST_LST_LSTCam", 1.0]], + "spe_width": [["type", "LST_LST_LSTCam", 1.0]], } } ) From 22e13001ea4f8a5290745e02232fcef96f95d723 Mon Sep 17 00:00:00 2001 From: Dan Parsons Date: Fri, 23 Jan 2026 14:58:17 +0100 Subject: [PATCH 23/43] Reverted change in iminuit requirement --- pyproject.toml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 5ed3b594b5b..b0229e991e4 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -57,7 +57,8 @@ all = [ "bokeh ~=3.0", "ctapipe[eventio]", "matplotlib ~=3.0", - "pyirf ~=0.13.0" + "pyirf ~=0.13.0", + "iminuit >=2" ] tests = [ @@ -100,8 +101,7 @@ dev = [ "setuptools_scm[toml]", ] -impact = ["tensorflow >=2.16.0", - "iminuit >=2"] +impact = ["tensorflow >=2.16.0"] [project.scripts] ctapipe-info = "ctapipe.tools.info:main" From fde474c6eae379f66ac3fb52b98b897b27127518 Mon Sep 17 00:00:00 2001 From: Dan Parsons Date: Fri, 23 Jan 2026 15:27:26 +0100 Subject: [PATCH 24/43] Fixin TF exceptions --- src/ctapipe/utils/template_network_interpolator.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/ctapipe/utils/template_network_interpolator.py b/src/ctapipe/utils/template_network_interpolator.py index 43eaedd0946..defd447eeb2 100644 --- a/src/ctapipe/utils/template_network_interpolator.py +++ b/src/ctapipe/utils/template_network_interpolator.py @@ -12,7 +12,8 @@ try: import tensorflow as tf except ModuleNotFoundError: - raise OptionalDependencyMissing("tensorflow") + tf = None + raise OptionalDependencyMissing("tensorflow") from None from .unstructured_interpolator import UnstructuredInterpolator From 46c97aa2b1e6e2560b052050113d713c9e935112 Mon Sep 17 00:00:00 2001 From: Dan Parsons Date: Tue, 27 Jan 2026 10:46:58 +0100 Subject: [PATCH 25/43] Added FreePACT tests and minor fixes to impact code. Added fallback to equivelent focal length if effective focal length is NaN. Removed mask on pixels below 0 intensity and replaced with a isclose to 0 statement. --- src/ctapipe/reco/freepact.py | 36 +++- src/ctapipe/reco/impact.py | 12 +- src/ctapipe/reco/tests/test_freepact.py | 194 ++++++++++++++++++ .../reco/tests/test_reconstruction_methods.py | 10 +- 4 files changed, 244 insertions(+), 8 deletions(-) create mode 100644 src/ctapipe/reco/tests/test_freepact.py diff --git a/src/ctapipe/reco/freepact.py b/src/ctapipe/reco/freepact.py index d40f7344064..fd0a5a7ec87 100644 --- a/src/ctapipe/reco/freepact.py +++ b/src/ctapipe/reco/freepact.py @@ -4,6 +4,7 @@ import numpy as np import numpy.ma as ma +import tensorflow as tf from ctapipe.core import traits from ctapipe.core.telescope_component import TelescopeParameter @@ -13,7 +14,7 @@ from ..compat import COPY_IF_NEEDED from .impact_utilities import rotate_translate -__all__ = ["FreePACTReconstructor"] +__all__ = ["FreePACTReconstructor", "create_dummy_freepact_templates"] class FreePACTReconstructor(ImPACTReconstructor): @@ -115,7 +116,6 @@ def get_likelihood( ) pix_x_rot = np.ma.array(pix_x_rot, mask=self.pixel_x.mask, copy=COPY_IF_NEEDED) pix_y_rot = np.ma.array(pix_y_rot, mask=self.pixel_y.mask, copy=COPY_IF_NEEDED) - # In the interpolator class we can gain speed advantages by using masked arrays # so we need to make sure here everything is masked likelihood = ma.zeros(self.image.shape) @@ -134,13 +134,14 @@ def get_likelihood( np.rad2deg(pix_y_rot[template_mask]), self.image[template_mask], ) + likelihood.mask = ma.getmask(self.image) + if goodness_of_fit: # return -2 * np.sum(likelihood[self.image>5]) / np.sum(self.image>5) # print(ma.getmask(self.image)) return -2 * ma.sum(likelihood) / np.sum(~ma.getmask(self.image)) likelihood = ma.sum(likelihood) * -2 - return likelihood @@ -160,3 +161,32 @@ class FreePACTProtonReconstructor(FreePACTReconstructor): def __init__(self, *args, **kwargs): super().__init__(*args, **kwargs) + + +def create_dummy_freepact_templates( + output_dir, telescope_type, zenith, azimuth, offset +): + """Create file with dummy template library + + Args: + output_dir (str): Output directory + telescope_type (str): Telescope type + zenith (float): Zenith angle in radians + azimuth (float): Azimuth angle in radians + offset (float): Offset in degrees + """ + model = tf.keras.Sequential( + [ + tf.keras.layers.InputLayer(shape=(6,)), + tf.keras.layers.Dense(5, activation="sigmoid", kernel_initializer="zeros"), + tf.keras.layers.Dense(5, activation="sigmoid", kernel_initializer="zeros"), + tf.keras.layers.Dense(1, activation="sigmoid", kernel_initializer="zeros"), + ] + ) + + model.compile(optimizer="adam", loss="mse") + + # reformat line to create file name string + output_file = f"/predict_{telescope_type}_{int(zenith)}deg_{int(azimuth)}deg_{offset:.1f}off.keras" + + model.save(output_dir + output_file) diff --git a/src/ctapipe/reco/impact.py b/src/ctapipe/reco/impact.py index 0c6641b1132..6ae9be9c3e9 100644 --- a/src/ctapipe/reco/impact.py +++ b/src/ctapipe/reco/impact.py @@ -590,6 +590,10 @@ class before minimisation can take place. This simply copies a mask = mask_dict[tel_id] focal_length = subarray.tel[tel_id].optics.effective_focal_length + if np.isnan(focal_length): + # fallback to equivalent focal length, do we warn here? + focal_length = subarray.tel[tel_id].optics.equivalent_focal_length + camera_frame = CameraFrame( telescope_pointing=telescope_pointing[tel_id], focal_length=focal_length, @@ -642,7 +646,8 @@ class before minimisation can take place. This simply copies a self.spe[i][:array_len] = p_spe[i] # Set the image mask - mask = self.image <= 0.0 + mask = np.isclose(self.image, 0, atol=1e-10, equal_nan=False) + self.pixel_x[mask], self.pixel_y[mask] = ma.masked, ma.masked self.image[mask] = ma.masked self.time[mask] = ma.masked @@ -736,7 +741,7 @@ def predict( if energy_reco == "FreePACTReconstructor": continue - energy = energy_seed[energy_reco].energy.value + energy = energy_seed[energy_reco].energy.to(u.TeV).value # for i in range(1): # energy = 0.5 # TEST!!! @@ -751,13 +756,12 @@ def predict( limits=limits, ) # fit_params_min, like = self.energy_guess(seed) - if like < like_min: + if like <= like_min: fit_params = fit_params_min like_min = like if fit_params is None: return self.INVALID_GEOMETRY, self.INVALID_ENERGY - # Now do full minimisation seed = create_seed( fit_params[0], fit_params[1], fit_params[2], fit_params[3], fit_params[4] diff --git a/src/ctapipe/reco/tests/test_freepact.py b/src/ctapipe/reco/tests/test_freepact.py new file mode 100644 index 00000000000..2c606c58a72 --- /dev/null +++ b/src/ctapipe/reco/tests/test_freepact.py @@ -0,0 +1,194 @@ +import os + +import astropy.units as u +import numpy as np +import pytest +from astropy.coordinates import AltAz, Angle +from traitlets.config import Config + +from ctapipe.containers import ( + HillasParametersContainer, + ReconstructedEnergyContainer, + ReconstructedGeometryContainer, +) +from ctapipe.reco.freepact import FreePACTReconstructor, create_dummy_freepact_templates +from ctapipe.utils import get_dataset_path + +pytest.importorskip("iminuit") +pytest.importorskip("tensorflow") + +SIMTEL_PATH = get_dataset_path( + "gamma_20deg_0deg_run2___cta-prod5-paranal_desert" + "-2147m-Paranal-dark_cone10-100evts.simtel.zst" +) + + +def get_simtel_profile_from_eventsource(): + """get a TableAtmosphereDensityProfile from a simtel file""" + from ctapipe.io import EventSource + + with EventSource(SIMTEL_PATH) as source: + return source.atmosphere_density_profile + + +@pytest.fixture(scope="session") +def table_profile(): + """a table profile for testing""" + return get_simtel_profile_from_eventsource() + + +class TestFreePACT: + @classmethod + def setup_class(self): + self.horizon_frame = AltAz() + + self.h1 = HillasParametersContainer( + fov_lon=1 * u.deg, + fov_lat=1 * u.deg, + r=1 * u.deg, + phi=Angle(0 * u.deg), + intensity=100, + length=0.4 * u.deg, + width=0.4 * u.deg, + psi=Angle(0 * u.deg), + skewness=0, + kurtosis=0, + ) + + def test_interpolation(self, tmp_path, example_subarray, table_profile): + """Test interpolation works on dummy template library""" + + for tel_type in example_subarray.telescope_types: + os.mkdir(str(tmp_path) + "/" + str(tel_type)[0:3] + "/") + create_dummy_freepact_templates( + str(tmp_path) + "/" + str(tel_type)[0:3] + "/", + str(tel_type)[0:3], + 0, + 0, + 0, + ) + create_dummy_freepact_templates( + str(tmp_path) + "/" + str(tel_type)[0:3] + "/", + str(tel_type)[0:3], + 20, + 0, + 0, + ) + create_dummy_freepact_templates( + str(tmp_path) + "/" + str(tel_type)[0:3] + "/", + str(tel_type)[0:3], + 0, + 180, + 0, + ) + create_dummy_freepact_templates( + str(tmp_path) + "/" + str(tel_type)[0:3] + "/", + str(tel_type)[0:3], + 20, + 180, + 0, + ) + + freepact_config = Config( + { + "FreePACTReconstructor": { + "image_template_path": [ + [ + "type", + "*", + str(tmp_path) + "/SST/", + ], + ] + } + } + ) + + freepact = FreePACTReconstructor( + example_subarray, table_profile, config=freepact_config + ) + + pred = freepact.prediction[list(freepact.prediction.keys())[0]]( + 0, + 0, + np.array([1, 1]), + np.array([100, 100]), + np.array([300, 300]), + np.array([[0, 0], [0, 0]]), + np.array([[0, 1], [0, 0]]), + np.array([[1, 1], [0, 0]]), + ) + assert pred.shape == (2, 2) + + +def test_selected_subarray( + subarray_and_event_gamma_off_axis_500_gev, tmp_path, table_profile +): + subarray, event = subarray_and_event_gamma_off_axis_500_gev + + for tel_type in subarray.telescope_types: + os.mkdir(str(tmp_path) + "/" + str(tel_type)[0:3] + "/") + create_dummy_freepact_templates( + str(tmp_path) + "/" + str(tel_type)[0:3] + "/", str(tel_type)[0:3], 0, 0, 0 + ) + create_dummy_freepact_templates( + str(tmp_path) + "/" + str(tel_type)[0:3] + "/", str(tel_type)[0:3], 60, 0, 0 + ) + create_dummy_freepact_templates( + str(tmp_path) + "/" + str(tel_type)[0:3] + "/", + str(tel_type)[0:3], + 0, + 180, + 0, + ) + create_dummy_freepact_templates( + str(tmp_path) + "/" + str(tel_type)[0:3] + "/", + str(tel_type)[0:3], + 60, + 180, + 0, + ) + print(str(tmp_path) + "/" + str(tel_type)[0:3] + "/") + freepact_config = Config( + { + "FreePACTReconstructor": { + "image_template_path": [ + [ + "type", + "LST_LST_LSTCam", + str(tmp_path) + "/LST/", + ], + ], + "pedestal_width": [["type", "LST_LST_LSTCam", 1.0]], + "spe_width": [["type", "LST_LST_LSTCam", 1.0]], + } + } + ) + + shower_test = ReconstructedGeometryContainer() + energy_test = ReconstructedEnergyContainer() + shower_test.prefix = "test" + energy_test.prefix = "test_energy" + # Transform everything back to a useful system + shower_test.alt, shower_test.az = 70 * u.deg, 0 * u.deg + + shower_test.core_x = -10 * u.m + shower_test.core_y = -10 * u.m + shower_test.core_tilted_x = -10 * u.m + shower_test.core_tilted_y = -10 * u.m + + shower_test.is_valid = True + + energy_test.energy = 0.5 * u.TeV + energy_test.is_valid = True + + event.dl2.stereo.geometry["test"] = shower_test + event.dl2.stereo.energy["test_energy"] = energy_test + + freepact_reco = FreePACTReconstructor( + subarray, table_profile, config=freepact_config + ) + freepact_reco(event) # effective focal length returning NaN + print(event.dl2.stereo.geometry["FreePACTReconstructor"]) + + assert event.dl2.stereo.geometry["FreePACTReconstructor"].is_valid + assert event.dl2.stereo.energy["FreePACTReconstructor"].is_valid diff --git a/src/ctapipe/reco/tests/test_reconstruction_methods.py b/src/ctapipe/reco/tests/test_reconstruction_methods.py index 1f8c8cefef1..af81e6dea72 100644 --- a/src/ctapipe/reco/tests/test_reconstruction_methods.py +++ b/src/ctapipe/reco/tests/test_reconstruction_methods.py @@ -62,7 +62,15 @@ def test_reconstructors(cls, tmp_path): "MST_MST_NectarCam", str(tmp_path) + "/MST_MST_NectarCam.template.gz", ], - ] + ], + "pedestal_width": [ + ["type", "LST_LST_LSTCam", 1.0], + ["type", "MST_MST_NectarCam", 1.0], + ], + "spe_width": [ + ["type", "LST_LST_LSTCam", 1.0], + ["type", "MST_MST_NectarCam", 1.0], + ], } } ) From 53007eb6e8f20820de2f4eb840323d24d960cfec Mon Sep 17 00:00:00 2001 From: Dan Parsons Date: Tue, 27 Jan 2026 13:24:01 +0100 Subject: [PATCH 26/43] Fixing import error as minor doc changes --- pyproject.toml | 4 +- src/ctapipe/reco/freepact.py | 32 ++++-- src/ctapipe/reco/impact.py | 3 - .../utils/template_network_interpolator.py | 107 ++++-------------- 4 files changed, 50 insertions(+), 96 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 442132e4129..2b7a3709e16 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -101,7 +101,9 @@ dev = [ "setuptools_scm[toml]", ] -impact = ["tensorflow >=2.16.0"] +impact = [ + "iminuit >=2", + "tensorflow >=2.16.0"] [project.scripts] ctapipe-info = "ctapipe.tools.info:main" diff --git a/src/ctapipe/reco/freepact.py b/src/ctapipe/reco/freepact.py index fd0a5a7ec87..93a90e8c290 100644 --- a/src/ctapipe/reco/freepact.py +++ b/src/ctapipe/reco/freepact.py @@ -19,9 +19,29 @@ class FreePACTReconstructor(ImPACTReconstructor): """ - FreePACT Reconstructor This class implements the FreePACT reconstructor, which uses a neural network - to predict the likelihood of camera images based on shower parameters. + to predict the likelihood of camera images based on shower parameters. The + reconstructor is based on the work presented in schwefer24. + + Because this application is computationally intensive the usual + advice to use astropy units for all quantities is ignored (as + these slow down some computations), instead units within the class + are fixed: + + - Angular units in radians + - Distance units in metres + - Energy units in TeV + + Parameters + ---------- + subarray : ctapipe.instrument.SubarrayDescription + The telescope subarray to use for reconstruction + atmosphere_profile : ctapipe.atmosphere.AtmosphereDensityProfile + Density vs. altitude profile of the local atmosphere + + References + ---------- + .. [schwefer24] Schwefer, Parsons, & Hinton, Astroparticle Physics 163 (2024), 103008 """ image_template_path = TelescopeParameter( @@ -32,9 +52,6 @@ class FreePACTReconstructor(ImPACTReconstructor): help=("Path to the image templates to be used in the reconstruction"), ).tag(config=True) - def __init__(self, *args, **kwargs): - super().__init__(*args, **kwargs) - def set_up_templates(self): """Set up the templates for the FreePACT reconstructor.""" template_sort_dict = {} @@ -146,7 +163,7 @@ def get_likelihood( class FreePACTProtonReconstructor(FreePACTReconstructor): - """FreePACT Proton Reconstructor + """ This class implements the FreePACT reconstructor for proton showers. It uses a neural network to predict the likelihood of camera images based on shower parameters. """ @@ -159,9 +176,6 @@ class FreePACTProtonReconstructor(FreePACTReconstructor): help=("Path to the image templates to be used in the reconstruction"), ).tag(config=True) - def __init__(self, *args, **kwargs): - super().__init__(*args, **kwargs) - def create_dummy_freepact_templates( output_dir, telescope_type, zenith, azimuth, offset diff --git a/src/ctapipe/reco/impact.py b/src/ctapipe/reco/impact.py index 6ae9be9c3e9..f2b9ddc3428 100644 --- a/src/ctapipe/reco/impact.py +++ b/src/ctapipe/reco/impact.py @@ -70,9 +70,6 @@ class ImPACTReconstructor(HillasGeometryReconstructor): templates to perform a maximum likelihood fit for the shower axis, energy and height of maximum. - Besides the image information, there is also the option to use the time gradient of the pixels - across the image as additional information in the fit. This requires an additional set of templates - Because this application is computationally intensive the usual advice to use astropy units for all quantities is ignored (as these slow down some computations), instead units within the class diff --git a/src/ctapipe/utils/template_network_interpolator.py b/src/ctapipe/utils/template_network_interpolator.py index defd447eeb2..7da573e6aee 100644 --- a/src/ctapipe/utils/template_network_interpolator.py +++ b/src/ctapipe/utils/template_network_interpolator.py @@ -11,13 +11,18 @@ try: import tensorflow as tf + + tf.config.set_visible_devices([], "GPU") + tf_function = tf.function except ModuleNotFoundError: tf = None - raise OptionalDependencyMissing("tensorflow") from None -from .unstructured_interpolator import UnstructuredInterpolator + # dummy decorator + def tf_function(f): + return f + -tf.config.set_visible_devices([], "GPU") +from .unstructured_interpolator import UnstructuredInterpolator class BaseTemplate: @@ -322,82 +327,6 @@ def __eq__(self, other): return self.template_file == other.template_file -class TimeGradientInterpolator(BaseTemplate): - """ - Class for interpolating between the time gradient predictions - """ - - def __init__(self, template_file): - """ - Parameters - ---------- - template_file: str - Location of pickle file containing ImPACT NN templates - """ - - super().__init__() - - self.template_file = template_file - - input_dict = None - with gzip.open(template_file, "r") as file_list: - input_dict = pickle.load(file_list) - - # The dict has at its highest level two keys: data and tel_type. - # The first just contains the template dict, the second a string to validate the tel_type - # that the templates are made for. - - data_input_dict = input_dict["data"] - self.tel_type_string = input_dict["tel_type"] - - keys = np.array(list(data_input_dict.keys())) - values = np.array(list(data_input_dict.values()), dtype=np.float32) - self.no_zenaz = False - - # First check if we even have a zen and azimuth entry - if len(keys[0]) > 4: - # If we do then for the sake of speed lets - self._create_table_matrix(keys, values) - else: - # If not we work as before - # Currently impact is not set up for offset dependent templates. - # Therefore remove offset (last) dimension from interpolator - self.interpolator = UnstructuredInterpolator( - keys[:-1], values, remember_last=False - ) - self.no_zenaz = True - - def __call__(self, zenith, azimuth, energy, impact, xmax): - """ - Evaluate expected time gradient for a set of shower parameters and pixel positions - Parameters - ---------- - energy: array-like - Energy of interpolated template - impact: array-like - Impact distance of interpolated template - xmax: array-like - Depth of maximum of interpolated templates - Returns - ------- - ndarray: Time Gradient expectation and RMS values - """ - array = np.stack((energy, impact, xmax), axis=-1) - - if self.no_zenaz: - interpolated_value = self.interpolator(array, None) - else: - interpolated_value = self.perform_interpolation( - zenith, azimuth, array, None - ) - - return interpolated_value - - # Two interpolators are the same if they are generated from the same file - def __eq__(self, other): - return self.template_file == other.template_file - - def load_prediction_files_filtered(directory): """ Finds all files in the given directory matching the pattern: @@ -451,7 +380,7 @@ def custom_symlog(value, linear_threshold=10.0): return value -@tf.function +@tf_function def evaluate_model(model, parameters): """ Evaluate the model with the given parameters. TF decorator used to @@ -466,13 +395,13 @@ def evaluate_model(model, parameters): Returns ------- - tf.Tensor + tf. Tensor The output of the model """ return model(parameters) -@tf.function +@tf_function def evaluate_model_interpolate( model_ul, model_uu, @@ -524,7 +453,7 @@ def evaluate_model_interpolate( v_ul = model_ul(parameters) v_lu = model_lu(parameters) v_uu = model_uu(parameters) - # print(v_ll, zenith) + v_l = ((v_ul - v_ll) * (zenith - zenith_l) / (zenith_u - zenith_l)) + v_ll v_u = ((v_uu - v_lu) * (zenith - zenith_l) / (zenith_u - zenith_l)) + v_lu @@ -688,6 +617,7 @@ def _evaluate_interpolator( """ if self.interpolator[zenith_bin][azimuth_bin] is None: self._create_interpolator(zenith_bin, azimuth_bin) + return np.asarray( evaluate_model(self.interpolator[zenith_bin][azimuth_bin], points) ) @@ -764,3 +694,14 @@ def __call__(self, energy, impact, xmax): def reset(self): """Reset the interpolator.""" return True + + +class DummyFreePACTInterpolator: + """Dummy FreePACT interpolator for testing purposes.""" + + def __call__(self, zenith, azimuth, energy, impact, xmax, xb, yb, amplitude): + return np.ones_like(xb) + + def reset(self): + """Reset the interpolator.""" + return True From e8e162fde4ab00f93f0c0c9babb2f7203590bdbf Mon Sep 17 00:00:00 2001 From: Dan Parsons Date: Tue, 27 Jan 2026 13:29:22 +0100 Subject: [PATCH 27/43] Again fixing optional dependencies for tensorflow. --- src/ctapipe/reco/freepact.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/src/ctapipe/reco/freepact.py b/src/ctapipe/reco/freepact.py index 93a90e8c290..9cfc52d304a 100644 --- a/src/ctapipe/reco/freepact.py +++ b/src/ctapipe/reco/freepact.py @@ -4,7 +4,11 @@ import numpy as np import numpy.ma as ma -import tensorflow as tf + +try: + import tensorflow as tf +except ModuleNotFoundError: + tf = None from ctapipe.core import traits from ctapipe.core.telescope_component import TelescopeParameter @@ -12,6 +16,7 @@ from ctapipe.utils.template_network_interpolator import FreePACTInterpolator from ..compat import COPY_IF_NEEDED +from ..exceptions import OptionalDependencyMissing from .impact_utilities import rotate_translate __all__ = ["FreePACTReconstructor", "create_dummy_freepact_templates"] @@ -189,6 +194,10 @@ def create_dummy_freepact_templates( azimuth (float): Azimuth angle in radians offset (float): Offset in degrees """ + + if tf is None: + raise OptionalDependencyMissing("tensorflow") + model = tf.keras.Sequential( [ tf.keras.layers.InputLayer(shape=(6,)), From aa90a60a712fad34a8069aafa640ba593d56edd0 Mon Sep 17 00:00:00 2001 From: Dan Parsons Date: Tue, 27 Jan 2026 14:00:45 +0100 Subject: [PATCH 28/43] Remove masked array copy function no longer supported in numpy 2 --- src/ctapipe/reco/impact.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/ctapipe/reco/impact.py b/src/ctapipe/reco/impact.py index f2b9ddc3428..153190244a4 100644 --- a/src/ctapipe/reco/impact.py +++ b/src/ctapipe/reco/impact.py @@ -14,7 +14,6 @@ from ctapipe.core.telescope_component import TelescopeParameter from ctapipe.exceptions import OptionalDependencyMissing -from ..compat import COPY_IF_NEEDED from ..containers import ReconstructedEnergyContainer, ReconstructedGeometryContainer from ..coordinates import ( CameraFrame, @@ -454,8 +453,8 @@ def get_likelihood( pix_x_rot, pix_y_rot = rotate_translate( self.pixel_y.data, self.pixel_x.data, source_y, source_x, -phi ) - pix_x_rot = np.ma.array(pix_x_rot, mask=self.pixel_x.mask, copy=COPY_IF_NEEDED) - pix_y_rot = np.ma.array(pix_y_rot, mask=self.pixel_y.mask, copy=COPY_IF_NEEDED) + pix_x_rot = np.ma.array(pix_x_rot, mask=self.pixel_x.mask) + pix_y_rot = np.ma.array(pix_y_rot, mask=self.pixel_y.mask) # In the interpolator class we can gain speed advantages by using masked arrays # so we need to make sure here everything is masked From da3af8fe023e49037be3a55e68baf0f09fce89a6 Mon Sep 17 00:00:00 2001 From: Dan Parsons Date: Tue, 27 Jan 2026 14:35:03 +0100 Subject: [PATCH 29/43] Remove masked array copy function no longer supported --- src/ctapipe/reco/freepact.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/ctapipe/reco/freepact.py b/src/ctapipe/reco/freepact.py index 9cfc52d304a..ec409c50d06 100644 --- a/src/ctapipe/reco/freepact.py +++ b/src/ctapipe/reco/freepact.py @@ -15,7 +15,6 @@ from ctapipe.reco.impact import ImPACTReconstructor from ctapipe.utils.template_network_interpolator import FreePACTInterpolator -from ..compat import COPY_IF_NEEDED from ..exceptions import OptionalDependencyMissing from .impact_utilities import rotate_translate @@ -136,8 +135,8 @@ def get_likelihood( pix_x_rot, pix_y_rot = rotate_translate( self.pixel_y.data, self.pixel_x.data, source_y, source_x, -phi ) - pix_x_rot = np.ma.array(pix_x_rot, mask=self.pixel_x.mask, copy=COPY_IF_NEEDED) - pix_y_rot = np.ma.array(pix_y_rot, mask=self.pixel_y.mask, copy=COPY_IF_NEEDED) + pix_x_rot = np.ma.array(pix_x_rot, mask=self.pixel_x.mask) + pix_y_rot = np.ma.array(pix_y_rot, mask=self.pixel_y.mask) # In the interpolator class we can gain speed advantages by using masked arrays # so we need to make sure here everything is masked likelihood = ma.zeros(self.image.shape) From c7025737a0c812b2c52f0836a48561cb10491bfd Mon Sep 17 00:00:00 2001 From: Dan Parsons Date: Tue, 27 Jan 2026 21:28:59 +0100 Subject: [PATCH 30/43] Add filter for tensorflow warnings --- pyproject.toml | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/pyproject.toml b/pyproject.toml index 2b7a3709e16..8770997b12c 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -177,7 +177,11 @@ filterwarnings = [ "ignore:'endQuoteChar':DeprecationWarning:matplotlib", "ignore:'unquoteResults':DeprecationWarning:matplotlib", "ignore:'parseAll':DeprecationWarning:pyparsing", + "ignore:__array__ implementation doesn't accept a copy keyword:DeprecationWarning", + "ignore:Type .* uses PyType_Spec with a metaclass that has custom tp_new:DeprecationWarning", + "ignore:datetime.datetime.utcfromtimestamp\\(\\):DeprecationWarning", ] + norecursedirs = [ ".git", "_build", From a396795ecb8ba02666a98d335d0be531f51e59a4 Mon Sep 17 00:00:00 2001 From: Dan Parsons Date: Wed, 28 Jan 2026 08:53:54 +0100 Subject: [PATCH 31/43] Filter for pyparsing warnings. --- pyproject.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/pyproject.toml b/pyproject.toml index 8770997b12c..fd909d4a858 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -180,6 +180,7 @@ filterwarnings = [ "ignore:__array__ implementation doesn't accept a copy keyword:DeprecationWarning", "ignore:Type .* uses PyType_Spec with a metaclass that has custom tp_new:DeprecationWarning", "ignore:datetime.datetime.utcfromtimestamp\\(\\):DeprecationWarning", + "ignore:pyparsing.warnings.PyparsingDeprecationWarning", ] norecursedirs = [ From 31ceb96c64553c3f81e8c81b845e62ad37d480a8 Mon Sep 17 00:00:00 2001 From: Dan Parsons Date: Wed, 28 Jan 2026 09:43:24 +0100 Subject: [PATCH 32/43] fix pyparsing warning --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index fd909d4a858..5c14c7a0d16 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -180,7 +180,7 @@ filterwarnings = [ "ignore:__array__ implementation doesn't accept a copy keyword:DeprecationWarning", "ignore:Type .* uses PyType_Spec with a metaclass that has custom tp_new:DeprecationWarning", "ignore:datetime.datetime.utcfromtimestamp\\(\\):DeprecationWarning", - "ignore:pyparsing.warnings.PyparsingDeprecationWarning", + "ignore:'enablePackrat' deprecated - use 'enable_packrat':pyparsing.warnings.PyparsingDeprecationWarning", ] norecursedirs = [ From 1fb8b0ee5ba0aada8ed7c58c4ac4dba2803da40c Mon Sep 17 00:00:00 2001 From: Dan Parsons Date: Wed, 28 Jan 2026 09:44:09 +0100 Subject: [PATCH 33/43] Again parsin warning fix --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 5c14c7a0d16..d91ef339144 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -180,7 +180,7 @@ filterwarnings = [ "ignore:__array__ implementation doesn't accept a copy keyword:DeprecationWarning", "ignore:Type .* uses PyType_Spec with a metaclass that has custom tp_new:DeprecationWarning", "ignore:datetime.datetime.utcfromtimestamp\\(\\):DeprecationWarning", - "ignore:'enablePackrat' deprecated - use 'enable_packrat':pyparsing.warnings.PyparsingDeprecationWarning", + "ignore::pyparsing.warnings.PyparsingDeprecationWarning", ] norecursedirs = [ From 25c0b89824befc7ee08ea6431150dc9af4150c39 Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Wed, 28 Jan 2026 09:41:33 +0100 Subject: [PATCH 34/43] Fix filterwarnings --- pyproject.toml | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index d91ef339144..a828786777e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -58,7 +58,7 @@ all = [ "ctapipe[eventio]", "matplotlib ~=3.0", "pyirf ~=0.13.0", - "iminuit >=2" + "iminuit >=2", ] tests = [ @@ -103,7 +103,8 @@ dev = [ impact = [ "iminuit >=2", - "tensorflow >=2.16.0"] + "tensorflow >=2.16.0", +] [project.scripts] ctapipe-info = "ctapipe.tools.info:main" @@ -177,10 +178,12 @@ filterwarnings = [ "ignore:'endQuoteChar':DeprecationWarning:matplotlib", "ignore:'unquoteResults':DeprecationWarning:matplotlib", "ignore:'parseAll':DeprecationWarning:pyparsing", + "ignore:'enablePackrat':DeprecationWarning:matplotlib", + # tensorflow (FreePACT) related warnings "ignore:__array__ implementation doesn't accept a copy keyword:DeprecationWarning", + # protobuf related warnings (via tensorflow) "ignore:Type .* uses PyType_Spec with a metaclass that has custom tp_new:DeprecationWarning", "ignore:datetime.datetime.utcfromtimestamp\\(\\):DeprecationWarning", - "ignore::pyparsing.warnings.PyparsingDeprecationWarning", ] norecursedirs = [ From d47afd10f7b4dfad584a265c44d716db20dcb760 Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Wed, 28 Jan 2026 10:19:34 +0100 Subject: [PATCH 35/43] Remove impact extra from minimal-dependencies CI job --- .github/workflows/ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 1c0e555f886..f9b2c163a51 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -51,7 +51,7 @@ jobs: os: ubuntu-24.04 python-version: "3.12" install-method: pip - extras: tests,impact + extras: tests - name: Linux (3.12, pip, coverage) os: ubuntu-24.04 From 6683ef3fdabc6c7de5a3238ff54524b53030cf6b Mon Sep 17 00:00:00 2001 From: Dan Parsons Date: Wed, 28 Jan 2026 10:19:47 +0100 Subject: [PATCH 36/43] sonarqube error fixes --- src/ctapipe/reco/freepact.py | 2 -- src/ctapipe/reco/impact.py | 11 +---------- 2 files changed, 1 insertion(+), 12 deletions(-) diff --git a/src/ctapipe/reco/freepact.py b/src/ctapipe/reco/freepact.py index ec409c50d06..98f9b44af80 100644 --- a/src/ctapipe/reco/freepact.py +++ b/src/ctapipe/reco/freepact.py @@ -158,8 +158,6 @@ def get_likelihood( likelihood.mask = ma.getmask(self.image) if goodness_of_fit: - # return -2 * np.sum(likelihood[self.image>5]) / np.sum(self.image>5) - # print(ma.getmask(self.image)) return -2 * ma.sum(likelihood) / np.sum(~ma.getmask(self.image)) likelihood = ma.sum(likelihood) * -2 diff --git a/src/ctapipe/reco/impact.py b/src/ctapipe/reco/impact.py index 153190244a4..854b4f61050 100644 --- a/src/ctapipe/reco/impact.py +++ b/src/ctapipe/reco/impact.py @@ -758,16 +758,7 @@ def predict( if fit_params is None: return self.INVALID_GEOMETRY, self.INVALID_ENERGY - # Now do full minimisation - seed = create_seed( - fit_params[0], fit_params[1], fit_params[2], fit_params[3], fit_params[4] - ) - - # fit_params, errors, like = self.minimise( - # params=seed[0], - # step=seed[1], - # limits=seed[2], - # ) + # We may want to now run a full minimisation, depending on if we ran the previous minimisation to completion # Create a container class for reconstructed shower From 64d93ca8d142d7c23fb7e3512fb89a3d5e59f641 Mon Sep 17 00:00:00 2001 From: Dan Parsons Date: Wed, 28 Jan 2026 17:23:54 +0100 Subject: [PATCH 37/43] Fix sonarqube complaints --- src/ctapipe/utils/template_network_interpolator.py | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) diff --git a/src/ctapipe/utils/template_network_interpolator.py b/src/ctapipe/utils/template_network_interpolator.py index 7da573e6aee..08e82c3e392 100644 --- a/src/ctapipe/utils/template_network_interpolator.py +++ b/src/ctapipe/utils/template_network_interpolator.py @@ -370,9 +370,7 @@ def custom_symlog(value, linear_threshold=10.0): """ for i in range(value.shape[0]): for j in range(value.shape[1]): - if np.abs(value[i][j]) < linear_threshold: - value[i][j] = value[i][j] - else: + if np.abs(value[i][j]) > linear_threshold: value[i][j] = np.sign(value[i][j]) * ( np.log2(np.abs(value[i][j] / linear_threshold)) + linear_threshold ) @@ -481,7 +479,6 @@ def __init__(self, directory): raise OptionalDependencyMissing("tensorflow") data_input_dict = load_prediction_files_filtered(directory) - # self.tel_type_string = telescope_type keys = np.array(list(data_input_dict.keys())) values = list(data_input_dict.values()) @@ -520,7 +517,7 @@ def _create_interpolator(self, zenith_bin, azimuth_bin): # Select these values from our range of keys selection = np.logical_and(self.keys.T[0] == zenith, self.keys.T[1] == azimuth) - bin_sel = np.where(selection)[0][0] + bin_sel = np.nonzero(selection)[0][0] # Create interpolator using this selection # Currently freepact is not set up for offset dependent templates. @@ -530,7 +527,6 @@ def _create_interpolator(self, zenith_bin, azimuth_bin): # We can now remove these entries. self.keys = self.keys[np.invert(selection)] del self.values[bin_sel] - # self.values = self.values[np.invert(selection)] def perform_interpolation(self, zenith, azimuth, interpolation_array, points): """ @@ -622,8 +618,6 @@ def _evaluate_interpolator( evaluate_model(self.interpolator[zenith_bin][azimuth_bin], points) ) - # return self.interpolator[zenith_bin][azimuth_bin].predict(interpolation_array, verbose=0, batch_size=10000) - def __call__(self, zenith, azimuth, energy, impact, xmax, xb, yb, amplitude): """ Evaluate expected time gradient for a set of shower parameters and pixel positions From 7e57590604e95b3d78ed58430d4b535f84f8cb2c Mon Sep 17 00:00:00 2001 From: Dan Parsons Date: Fri, 30 Jan 2026 14:19:07 +0100 Subject: [PATCH 38/43] First FreePACT trainer commit --- src/ctapipe/tools/train_freepact.py | 432 ++++++++++++++++++ .../utils/template_network_interpolator.py | 3 +- 2 files changed, 434 insertions(+), 1 deletion(-) create mode 100644 src/ctapipe/tools/train_freepact.py diff --git a/src/ctapipe/tools/train_freepact.py b/src/ctapipe/tools/train_freepact.py new file mode 100644 index 00000000000..6ef98df6f64 --- /dev/null +++ b/src/ctapipe/tools/train_freepact.py @@ -0,0 +1,432 @@ +""" +Tool for training a Pixel Classifier using a Keras MLP. +""" + +import astropy.units as u +import numpy as np +from tensorflow.keras.callbacks import EarlyStopping, ReduceLROnPlateau +from tensorflow.keras.layers import Dense +from tensorflow.keras.models import Sequential +from tensorflow.keras.optimizers import Adam + +from ctapipe.core import Tool +from ctapipe.core.traits import Int, IntTelescopeParameter, Path +from ctapipe.exceptions import InputMissing +from ctapipe.image import dilate, tailcuts_clean +from ctapipe.io import TableLoader +from ctapipe.utils.template_network_interpolator import custom_symlog + +__all__ = ["TrainFreePACT"] + + +class TrainFreePACT(Tool): + """ + Tool to train a pixel classifier (FreePACT) MLP on dl1 images. + + The tool loads DL1 images and true shower parameters, constructs a + pixel-wise dataset (Signal vs Background), and trains an MLP classifier. + Background is generated by shuffling amplitudes. + """ + + name = "ctapipe-train-freepact" + description = __doc__ + + output_path = Path( + directory_ok=False, + help="Output path for the trained model (Keras format).", + ).tag(config=True) + + n_events = IntTelescopeParameter( + default_value=None, + allow_none=True, + help=( + "Number of events to use for training per telescope type." + " If not given, all available events will be used." + ), + ).tag(config=True) + + batch_size = Int( + default_value=100000, + help="Batch size for training.", + ).tag(config=True) + + epochs = Int( + default_value=100, + help="Number of training epochs.", + ).tag(config=True) + + chunk_size = Int( + default_value=1000, + help="Number of events to load at once.", + ).tag(config=True) + + limit = Int( + default_value=None, + allow_none=True, + help="Limit on images used in training", + ).tag(config=True) + + aliases = { + ("i", "input"): "TableLoader.input_url", + ("o", "output"): "TrainFreePACT.output_path", + "epochs": "TrainFreePACT.epochs", + "n_events": "TrainFreePACT.n_events", + "batch_size": "TrainFreePACT.batch_size", + "chunk_size": "TrainFreePACT.chunk_size", + "imagelimit": "TrainFreePACT.limit", + } + + classes = [TableLoader] + + def setup(self): + """Initialize components.""" + # Check for tensorflow + try: + import tensorflow # noqa: F401 + except ImportError: + self.log.critical("tensorflow is required for this tool.") + self.exit(1) + + try: + self.loader = self.enter_context( + TableLoader( + parent=self, + dl1_images=True, + simulated=True, + true_parameters=True, + instrument=True, + ) + ) + except InputMissing: + self.log.critical("Specifying TableLoader.input_url is required.") + self.exit(1) + + if self.output_path is None: + self.log.critical("Output path is required via -o / --output.") + self.exit(1) + + def start(self): + """Load data and train models.""" + # We will train one model per telescope type + + types = self.loader.subarray.telescope_types + self.log.info("Found telescope types: %s", types) + + for tel_type in types: + self.log.info("Processing telescope type: %s", tel_type) + X, y = self._load_and_process_data(tel_type) + + if len(X) == 0: + self.log.warning("No data found for %s", tel_type) + continue + + self.log.info( + "Training model for %s with %d samples (batch_size=%d, epochs=%d)", + tel_type, + len(X), + self.batch_size, + self.epochs, + ) + + model = self._build_model(input_shape=(X.shape[1],)) + callbacks = [ + EarlyStopping( + monitor="val_loss", patience=50, verbose=1, start_from_epoch=50 + ), + ReduceLROnPlateau( + monitor="val_loss", + factor=0.5, + patience=20, + verbose=1, + start_from_epoch=50, + ), + ] + model.fit( + X, + y, + epochs=self.epochs, + batch_size=self.batch_size, + validation_split=0.1, + verbose=1, + callbacks=callbacks, + ) + + # Save model + # Construct a filename for this telescope type + # If the user provided a file path (e.g. model.keras), we might need to split it + # or save a dictionary of models? + # Keras models are usually saved as directories or files. + # Let's append the telescope type to the output name. + out_name = str(self.output_path).replace(".keras", "").replace(".h5", "") + save_path = f"{out_name}_{tel_type}.keras" + self.log.info("Saving model to %s", save_path) + model.save(save_path) + + def _load_and_process_data(self, tel_type): + """Read chunked data and create pixel-wise dataset.""" + from astropy.coordinates import AltAz, SkyCoord + + from ctapipe.coordinates import GroundFrame, NominalFrame, TiltedGroundFrame + + # Get camera geometry to map pixels to positions + example_tel_id = None + for tel_id, tel in self.loader.subarray.tel.items(): + if str(tel) == str(tel_type): + example_tel_id = tel_id + break + + if example_tel_id is None: + # Try matching by string representation if direct comparison failed + for tel_id, tel in self.loader.subarray.tel.items(): + if str(tel) == str(tel_type): + example_tel_id = tel_id + break + + if example_tel_id is None: + self.log.error("Could not find telescope for type %s", tel_type) + return np.array([]), np.array([]) + + geometry = self.loader.subarray.tel[example_tel_id].camera.geometry + focal_length = self.loader.subarray.tel[ + example_tel_id + ].optics.equivalent_focal_length + + pix_x = geometry.pix_x + pix_y = geometry.pix_y + + # Convert pixels to Nominal Frame (Field of View) + # Using a simple approximation where Nominal Frame is centered on the telescope optical axis + # and aligned with the camera. This neglects mispointing of the telescope relative to the array. + # If full transformation is needed, we would need to rotate for each event. + # For training a pixel classifier on single telescope images, + # using the telescope field of view coordinates is usually what is meant by "nominal frame". + + fov_lon = np.arctan(pix_x / focal_length).to_value(u.deg) + fov_lat = np.arctan(pix_y / focal_length).to_value(u.deg) + + # Iterate over chunks + iterator = self.loader.read_telescope_events_chunked( + chunk_size=self.chunk_size, + telescopes=[tel_type], + dl1_images=True, + dl1_parameters=False, + simulated=True, + true_parameters=True, + instrument=True, + pointing=True, + ) + + pixel_features_list = [] + n_loaded = 0 + + for chunk in iterator: + data = chunk.data + if len(data) == 0: + continue + + # Extract features + images = data["image"] + mask_array = np.zeros(images.shape, dtype=bool) + + for i in range(len(images)): + mask = tailcuts_clean(geometry, images[i], 10, 5) + if images[i].sum() < 60: + mask_array[i] = False + continue + for j in range(2): + mask = dilate(geometry, mask) + mask_array[i] = mask + + energies = data["true_energy"].quantity.to_value(u.TeV) + core_x = data["true_core_x"].quantity + core_y = data["true_core_y"].quantity + + # X max + if "true_x_max" in data.colnames: + x_max = data["true_x_max"].quantity.to_value(u.g / (u.cm**2)) + else: + x_max = np.zeros(len(data)) + if "true_h_max" not in data.colnames: # Avoid redundant warnings + self.log.warning("true_x_max not found, using 0") + + # Telescope position + # TableLoader with instrument=True usually adds pos_x, pos_y, pos_z in GroundFrame + if "pos_x" in data.colnames: + tel_x = data["pos_x"].quantity + tel_y = data["pos_y"].quantity + tel_z = data["pos_z"].quantity + else: + # Fallback to older names if necessary, but user suggested 'pos_x' + tel_x = data["tel_pos_x"].quantity + tel_y = data["tel_pos_y"].quantity + tel_z = data["tel_pos_z"].quantity + + # Tilted Frame Transformation + # correct core and telescope positions to be in the tilted frame + pointing_alt = data["telescope_pointing_altitude"].quantity + pointing_az = data["telescope_pointing_azimuth"].quantity + + # Create SkyCoord for pointing (one per event) + # TiltedGroundFrame needs a pointing direction. + + p_dir = AltAz(alt=pointing_alt, az=pointing_az) + tilted_frame = TiltedGroundFrame(pointing_direction=p_dir) + nominal_frame = NominalFrame(origin=p_dir) + + source_alt = data["true_alt"].quantity + source_az = data["true_az"].quantity + source_sky_coord = SkyCoord(alt=source_alt, az=source_az, frame=AltAz()) + source_nominal = source_sky_coord.transform_to(nominal_frame) + + source_x = source_nominal.fov_lon.to_value(u.radian) + source_y = source_nominal.fov_lat.to_value(u.radian) + + # Transform Core + # GroundFrame -> TiltedGroundFrame + # We assume core_z = 0 for the core position on ground? + # Usually true_core is on the ground plane (z=0 in GroundFrame). + core_ground = SkyCoord(x=core_x, y=core_y, z=0 * u.m, frame=GroundFrame()) + core_tilted = core_ground.transform_to(tilted_frame) + core_acc_x = core_tilted.x.to_value(u.m) + core_acc_y = core_tilted.y.to_value(u.m) + + # Transform Telescope Position + tel_ground = SkyCoord(x=tel_x, y=tel_y, z=tel_z, frame=GroundFrame()) + tel_tilted = tel_ground.transform_to(tilted_frame) + tel_acc_x = tel_tilted.x.to_value(u.m) + tel_acc_y = tel_tilted.y.to_value(u.m) + + # Calculate impact distance in Tilted Frame + dx = core_acc_x - tel_acc_x + dy = core_acc_y - tel_acc_y + phi = np.arctan2(dy, dx) + impact = np.sqrt(dx**2 + dy**2) + + # Memory Optimization: Use mask indices instead of tiling + # Get indices of valid pixels (rows=event_idx, cols=pixel_idx) + rows, cols = np.nonzero(mask_array) + + # Select features for valid pixels + energies_sel = energies[rows] + impact_sel = impact[rows] + xmax_sel = x_max[rows] + + # Select coordinates (Nominal Frame / FOV) + # fov_lon/lat are 1D arrays of pixel positions (Radians) + pix_lon_sel = fov_lon[cols] + pix_lat_sel = fov_lat[cols] + + # Select Source Position and Rotation Angle + # source_x/y and phi are per event (Radians) + source_x_sel = source_x[rows] + source_y_sel = source_y[rows] + phi_sel = phi[rows] + + # Rotation center on source + diff_lon = pix_lon_sel - source_x_sel + diff_lat = pix_lat_sel - source_y_sel + + # Rotate + cos_phi = np.cos(-phi_sel) + sin_phi = np.sin(-phi_sel) + + pix_lon_rot = diff_lon * cos_phi - diff_lat * sin_phi + pix_lat_rot = diff_lon * sin_phi + diff_lat * cos_phi + + # Select Images (Amplitudes) + images_sel = images[mask_array] + + # custom_symlog(np.array([images_sel/ 10]).ravel(), linear_threshold=2) + # tf.config.set_visible_devices([], "GPU") + # Stack features + X_chunk = np.stack( + [ + pix_lon_rot, + pix_lat_rot, + np.log10(energies_sel), + impact_sel / 100, + xmax_sel / 100, + custom_symlog( + np.array([images_sel / 10]), linear_threshold=2 + ).ravel(), + ], + axis=1, + ).astype(np.float32) + + pixel_features_list.append(X_chunk) + n_ev = np.sum(mask_array) + + n_loaded += n_ev + if self.limit and n_loaded >= self.limit: + break + + if not pixel_features_list: + return np.array([]), np.array([]) + + X_all = np.concatenate(pixel_features_list) + print(X_all.shape) + # Dataset creation (Signal vs Background) + n_total_pixels = len(X_all) + n_samples = n_total_pixels // 2 + + indices = np.random.choice(n_total_pixels, size=n_total_pixels, replace=False) + indices_sig = indices[:n_samples] + indices_bg = indices[n_samples : 2 * n_samples] + + X_sig = X_all[indices_sig] + y_sig = np.ones(len(X_sig)) + + X_bg = X_all[indices_bg].copy() + np.random.shuffle(X_bg[:, 5]) # Shuffle amplitude + y_bg = np.zeros(len(X_bg)) + + X = np.concatenate([X_sig, X_bg]) + y = np.concatenate([y_sig, y_bg]) + + train_idx = np.random.permutation(len(X)) + X = X[train_idx] + y = y[train_idx] + + return X, y + + def _build_model(self, input_shape): + """Build Keras MLP.""" + from tensorflow.keras.layers import Input + + model = Sequential( + [ + Input(shape=input_shape), + Dense(64, activation="swish"), + Dense(64, activation="swish"), + Dense(64, activation="swish"), + Dense(64, activation="swish"), + Dense(64, activation="swish"), + Dense(64, activation="swish"), + Dense(64, activation="swish"), + Dense(64, activation="swish"), + Dense(64, activation="swish"), + Dense(64, activation="swish"), + Dense(64, activation="swish"), + Dense(1, activation="sigmoid"), + ] + ) + + model.compile( + optimizer=Adam(learning_rate=0.001), + loss="binary_crossentropy", + metrics=["accuracy"], + ) + return model + + def finish(self): + """Cleanup.""" + self.loader.close() + + +def main(): + tool = TrainFreePACT() + tool.run() + + +if __name__ == "__main__": + main() diff --git a/src/ctapipe/utils/template_network_interpolator.py b/src/ctapipe/utils/template_network_interpolator.py index 08e82c3e392..976149c72b1 100644 --- a/src/ctapipe/utils/template_network_interpolator.py +++ b/src/ctapipe/utils/template_network_interpolator.py @@ -12,7 +12,6 @@ try: import tensorflow as tf - tf.config.set_visible_devices([], "GPU") tf_function = tf.function except ModuleNotFoundError: tf = None @@ -477,6 +476,8 @@ def __init__(self, directory): super().__init__() if tf is None: raise OptionalDependencyMissing("tensorflow") + else: + tf.config.set_visible_devices([], "GPU") data_input_dict = load_prediction_files_filtered(directory) From 85a2b966ad1b25c09a36f2563f8e2971d19ebcc5 Mon Sep 17 00:00:00 2001 From: Dan Parsons Date: Thu, 5 Feb 2026 10:29:58 +0100 Subject: [PATCH 39/43] Reordered parameters in the interpolator --- src/ctapipe/utils/template_network_interpolator.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ctapipe/utils/template_network_interpolator.py b/src/ctapipe/utils/template_network_interpolator.py index 976149c72b1..1f352d3e3d6 100644 --- a/src/ctapipe/utils/template_network_interpolator.py +++ b/src/ctapipe/utils/template_network_interpolator.py @@ -646,8 +646,8 @@ def __call__(self, zenith, azimuth, energy, impact, xmax, xb, yb, amplitude): ( xb.ravel(), yb.ravel(), - np.repeat(impact, repeat_num, axis=-1), np.repeat(energy, repeat_num, axis=-1), + np.repeat(impact, repeat_num, axis=-1), np.repeat(xmax, repeat_num, axis=-1), amplitude.ravel(), ), From 5b30303763e6aa5dddd82a0dcca8489f93a15c18 Mon Sep 17 00:00:00 2001 From: Dan Parsons Date: Thu, 5 Feb 2026 15:37:02 +0100 Subject: [PATCH 40/43] General refactor of training script to work with DL1 exported data --- src/ctapipe/tools/train_freepact.py | 507 +++++++++++++++++----------- 1 file changed, 301 insertions(+), 206 deletions(-) diff --git a/src/ctapipe/tools/train_freepact.py b/src/ctapipe/tools/train_freepact.py index 6ef98df6f64..1f70db9484e 100644 --- a/src/ctapipe/tools/train_freepact.py +++ b/src/ctapipe/tools/train_freepact.py @@ -4,14 +4,13 @@ import astropy.units as u import numpy as np -from tensorflow.keras.callbacks import EarlyStopping, ReduceLROnPlateau +from tensorflow.keras.callbacks import EarlyStopping, ModelCheckpoint, ReduceLROnPlateau from tensorflow.keras.layers import Dense from tensorflow.keras.models import Sequential from tensorflow.keras.optimizers import Adam from ctapipe.core import Tool -from ctapipe.core.traits import Int, IntTelescopeParameter, Path -from ctapipe.exceptions import InputMissing +from ctapipe.core.traits import Float, Int, IntTelescopeParameter, List, Path from ctapipe.image import dilate, tailcuts_clean from ctapipe.io import TableLoader from ctapipe.utils.template_network_interpolator import custom_symlog @@ -56,7 +55,7 @@ class TrainFreePACT(Tool): ).tag(config=True) chunk_size = Int( - default_value=1000, + default_value=10000, help="Number of events to load at once.", ).tag(config=True) @@ -66,14 +65,61 @@ class TrainFreePACT(Tool): help="Limit on images used in training", ).tag(config=True) + reweight_index = Float( + default_value=0, + allow_none=True, + help="Index for reweighting events by energy", + ).tag(config=True) + + reweight_energy = Float( + default_value=1, + allow_none=True, + help="Energy for reweighting events by energy", + ).tag(config=True) + + input_url = List( + Path(exists=True, directory_ok=False), + default_value=[], + help="Input DL1 files", + ).tag(config=True) + + layer_number = Int( + default_value=10, + help="Number of hidden layers in the MLP.", + ).tag(config=True) + + layer_neurons = Int( + default_value=50, + help="Number of neurons in each hidden layer.", + ).tag(config=True) + + tailcut_threshold1 = Float( + default_value=10, + help="Threshold 1 for tailcut cleaning.", + ).tag(config=True) + + tailcut_threshold2 = Float( + default_value=5, + help="Threshold 2 for tailcut cleaning.", + ).tag(config=True) + + tailcut_additional_rows = Int( + default_value=4, + help="Additional rows to clean around the core.", + ).tag(config=True) + aliases = { - ("i", "input"): "TableLoader.input_url", + ("i", "input"): "TrainFreePACT.input_url", ("o", "output"): "TrainFreePACT.output_path", "epochs": "TrainFreePACT.epochs", "n_events": "TrainFreePACT.n_events", "batch_size": "TrainFreePACT.batch_size", "chunk_size": "TrainFreePACT.chunk_size", "imagelimit": "TrainFreePACT.limit", + "reweight_index": "TrainFreePACT.reweight_index", + "reweight_energy": "TrainFreePACT.reweight_energy", + "layer_number": "TrainFreePACT.layer_number", + "layer_neurons": "TrainFreePACT.layer_neurons", } classes = [TableLoader] @@ -87,18 +133,24 @@ def setup(self): self.log.critical("tensorflow is required for this tool.") self.exit(1) + if not self.input_url: + self.log.critical("Specifying input file(s) is required.") + self.exit(1) + + # Open the first file to get the subarray description + # We assume all files have the same subarray try: - self.loader = self.enter_context( - TableLoader( - parent=self, - dl1_images=True, - simulated=True, - true_parameters=True, - instrument=True, - ) - ) - except InputMissing: - self.log.critical("Specifying TableLoader.input_url is required.") + with TableLoader( + input_url=self.input_url[0], + parent=self, + dl1_images=True, + simulated=True, + true_parameters=True, + instrument=True, + ) as loader: + self.subarray = loader.subarray + except Exception as e: + self.log.critical(f"Could not read subarray from first file: {e}") self.exit(1) if self.output_path is None: @@ -109,7 +161,7 @@ def start(self): """Load data and train models.""" # We will train one model per telescope type - types = self.loader.subarray.telescope_types + types = self.subarray.telescope_types self.log.info("Found telescope types: %s", types) for tel_type in types: @@ -128,7 +180,14 @@ def start(self): self.epochs, ) - model = self._build_model(input_shape=(X.shape[1],)) + model = self._build_model( + input_shape=(X.shape[1],), + layer_number=self.layer_number, + layer_neurons=self.layer_neurons, + ) + out_name = str(self.output_path).replace(".keras", "").replace(".h5", "") + save_path = f"{out_name}_{tel_type}.keras" + callbacks = [ EarlyStopping( monitor="val_loss", patience=50, verbose=1, start_from_epoch=50 @@ -140,6 +199,12 @@ def start(self): verbose=1, start_from_epoch=50, ), + ModelCheckpoint( + filepath=save_path, + monitor="val_loss", + save_best_only=True, + verbose=1, + ), ] model.fit( X, @@ -157,8 +222,7 @@ def start(self): # or save a dictionary of models? # Keras models are usually saved as directories or files. # Let's append the telescope type to the output name. - out_name = str(self.output_path).replace(".keras", "").replace(".h5", "") - save_path = f"{out_name}_{tel_type}.keras" + self.log.info("Saving model to %s", save_path) model.save(save_path) @@ -166,18 +230,23 @@ def _load_and_process_data(self, tel_type): """Read chunked data and create pixel-wise dataset.""" from astropy.coordinates import AltAz, SkyCoord - from ctapipe.coordinates import GroundFrame, NominalFrame, TiltedGroundFrame + from ctapipe.coordinates import ( + CameraFrame, + GroundFrame, + NominalFrame, + TiltedGroundFrame, + ) # Get camera geometry to map pixels to positions example_tel_id = None - for tel_id, tel in self.loader.subarray.tel.items(): + for tel_id, tel in self.subarray.tel.items(): if str(tel) == str(tel_type): example_tel_id = tel_id break if example_tel_id is None: # Try matching by string representation if direct comparison failed - for tel_id, tel in self.loader.subarray.tel.items(): + for tel_id, tel in self.subarray.tel.items(): if str(tel) == str(tel_type): example_tel_id = tel_id break @@ -186,177 +255,213 @@ def _load_and_process_data(self, tel_type): self.log.error("Could not find telescope for type %s", tel_type) return np.array([]), np.array([]) - geometry = self.loader.subarray.tel[example_tel_id].camera.geometry - focal_length = self.loader.subarray.tel[ - example_tel_id - ].optics.equivalent_focal_length + geometry = self.subarray.tel[example_tel_id].camera.geometry + focal_length = self.subarray.tel[example_tel_id].optics.effective_focal_length pix_x = geometry.pix_x pix_y = geometry.pix_y - # Convert pixels to Nominal Frame (Field of View) - # Using a simple approximation where Nominal Frame is centered on the telescope optical axis - # and aligned with the camera. This neglects mispointing of the telescope relative to the array. - # If full transformation is needed, we would need to rotate for each event. - # For training a pixel classifier on single telescope images, - # using the telescope field of view coordinates is usually what is meant by "nominal frame". - - fov_lon = np.arctan(pix_x / focal_length).to_value(u.deg) - fov_lat = np.arctan(pix_y / focal_length).to_value(u.deg) - - # Iterate over chunks - iterator = self.loader.read_telescope_events_chunked( - chunk_size=self.chunk_size, - telescopes=[tel_type], - dl1_images=True, - dl1_parameters=False, - simulated=True, - true_parameters=True, - instrument=True, - pointing=True, - ) - pixel_features_list = [] n_loaded = 0 - for chunk in iterator: - data = chunk.data - if len(data) == 0: - continue + for input_file in self.input_url: + self.log.info("Loading data from %s", input_file) + + with TableLoader( + input_url=input_file, + parent=self, + dl1_images=True, + simulated=True, + true_parameters=True, + instrument=True, + pointing=True, + dl1_parameters=False, + ) as loader: + # Iterate over chunks + iterator = loader.read_telescope_events_chunked( + chunk_size=self.chunk_size, + telescopes=[tel_type], + dl1_images=True, + dl1_parameters=False, + simulated=True, + true_parameters=True, + instrument=True, + pointing=True, + ) + + for chunk in iterator: + data = chunk.data + if len(data) == 0: + continue + + # Extract features + images = data["image"] + mask_array = np.zeros(images.shape, dtype=bool) + + for i in range(len(images)): + mask = tailcuts_clean( + geometry, + images[i], + self.tailcut_threshold1, + self.tailcut_threshold2, + ) + for j in range(self.tailcut_additional_rows): + mask = dilate(geometry, mask) + mask_array[i] = mask + + energies = data["true_energy"].quantity.to_value(u.TeV) + if self.reweight_index != 0: + weight_factor = ( + energies / (self.reweight_energy * u.TeV) + ) ** self.reweight_index + random_value = np.random.rand(mask_array.shape[0]) + reweight_mask = random_value < weight_factor.value + mask_array = np.logical_and( + mask_array, reweight_mask[:, np.newaxis] + ) + + core_x = data["true_core_x"].quantity + core_y = data["true_core_y"].quantity + true_zenith = 90 * u.deg - data["true_alt"].quantity + + # X max + if "true_x_max" in data.colnames: + x_max = data["true_x_max"].quantity.to_value(u.g / (u.cm**2)) + x_max /= np.cos(true_zenith) # convert to slant depth + else: + x_max = np.zeros(len(data)) + if ( + "true_h_max" not in data.colnames + ): # Avoid redundant warnings + self.log.warning("true_x_max not found, using 0") + + # Telescope position + # TableLoader with instrument=True usually adds pos_x, pos_y, pos_z in GroundFrame + if "pos_x" in data.colnames: + tel_x = data["pos_x"].quantity + tel_y = data["pos_y"].quantity + tel_z = data["pos_z"].quantity + else: + # Fallback to older names if necessary, but user suggested 'pos_x' + tel_x = data["tel_pos_x"].quantity + tel_y = data["tel_pos_y"].quantity + tel_z = data["tel_pos_z"].quantity + + # Tilted Frame Transformation + # correct core and telescope positions to be in the tilted frame + pointing_alt = data["telescope_pointing_altitude"].quantity + pointing_az = data["telescope_pointing_azimuth"].quantity + + p_dir = AltAz(alt=pointing_alt, az=pointing_az) + tilted_frame = TiltedGroundFrame(pointing_direction=p_dir) + # nominal_frame not needed here specifically, we do it per pixel + + # Transform Core + # GroundFrame -> TiltedGroundFrame + # We assume core_z = 0 for the core position on ground? + # Usually true_core is on the ground plane (z=0 in GroundFrame). + core_ground = SkyCoord( + x=core_x, y=core_y, z=0 * u.m, frame=GroundFrame() + ) + core_tilted = core_ground.transform_to(tilted_frame) + core_acc_x = core_tilted.x.to_value(u.m) + core_acc_y = core_tilted.y.to_value(u.m) + + # Transform Telescope Position + tel_ground = SkyCoord( + x=tel_x, y=tel_y, z=tel_z, frame=GroundFrame() + ) + tel_tilted = tel_ground.transform_to(tilted_frame) + tel_acc_x = tel_tilted.x.to_value(u.m) + tel_acc_y = tel_tilted.y.to_value(u.m) + + # Calculate impact distance in Tilted Frame + dx = core_acc_x - tel_acc_x + dy = core_acc_y - tel_acc_y + phi = np.arctan2(dy, dx) + impact = np.sqrt(dx**2 + dy**2) + + # Memory Optimization: Use mask indices instead of tiling + # Get indices of valid pixels (rows=event_idx, cols=pixel_idx) + rows, cols = np.nonzero(mask_array) + + # Select features for valid pixels + energies_sel = energies[rows] + impact_sel = impact[rows] + xmax_sel = x_max[rows] + + # Transform Camera Coordinates to Nominal Frame for valid pixels + # We need to construct the pointing for each valid pixel + pointing_alt_sel = pointing_alt[rows] + pointing_az_sel = pointing_az[rows] + pointing_sel = AltAz(alt=pointing_alt_sel, az=pointing_az_sel) + + # Construct frames for each valid pixel (vectorized) + # Note: CameraFrame and NominalFrame can broadcast their parameters + cam_frame = CameraFrame( + focal_length=focal_length, telescope_pointing=pointing_sel + ) + nom_frame = NominalFrame(origin=pointing_sel) + + pixel_coords = SkyCoord( + x=pix_x[cols], y=pix_y[cols], frame=cam_frame + ) + pixel_nominal = pixel_coords.transform_to(nom_frame) + + pix_lon_sel = pixel_nominal.fov_lon.to_value(u.deg) + pix_lat_sel = pixel_nominal.fov_lat.to_value(u.deg) + + # Select Source Position and Rotation Angle + # source_x/y and phi are per event (Radians) + # We need source position in Nominal Frame too to calculate rotation + source_alt_sel = data["true_alt"].quantity[rows] + source_az_sel = data["true_az"].quantity[rows] + source_sky_coord_sel = SkyCoord( + alt=source_alt_sel, az=source_az_sel, frame=AltAz() + ) + source_nominal_sel = source_sky_coord_sel.transform_to(nom_frame) + + source_x_sel = source_nominal_sel.fov_lon.to_value(u.deg) + source_y_sel = source_nominal_sel.fov_lat.to_value(u.deg) + + phi_sel = phi[rows] + + # Rotation center on source + diff_lon = pix_lon_sel - source_x_sel + diff_lat = pix_lat_sel - source_y_sel + + # Rotate + cos_phi = np.cos(-phi_sel) + sin_phi = np.sin(-phi_sel) + + pix_lon_rot = diff_lat * cos_phi - diff_lon * sin_phi + pix_lat_rot = diff_lat * sin_phi + diff_lon * cos_phi + + # Select Images (Amplitudes) + images_sel = images[mask_array] + + # Stack features + X_chunk = np.stack( + [ + pix_lon_rot, + pix_lat_rot, + np.log10(energies_sel), + impact_sel / 100, + xmax_sel / 100, + custom_symlog( + np.array([images_sel / 10]), linear_threshold=2 + ).ravel(), + ], + axis=1, + ).astype(np.float32) + + pixel_features_list.append(X_chunk) + n_ev = np.sum(mask_array) + + n_loaded += n_ev + if self.limit and n_loaded >= self.limit: + break - # Extract features - images = data["image"] - mask_array = np.zeros(images.shape, dtype=bool) - - for i in range(len(images)): - mask = tailcuts_clean(geometry, images[i], 10, 5) - if images[i].sum() < 60: - mask_array[i] = False - continue - for j in range(2): - mask = dilate(geometry, mask) - mask_array[i] = mask - - energies = data["true_energy"].quantity.to_value(u.TeV) - core_x = data["true_core_x"].quantity - core_y = data["true_core_y"].quantity - - # X max - if "true_x_max" in data.colnames: - x_max = data["true_x_max"].quantity.to_value(u.g / (u.cm**2)) - else: - x_max = np.zeros(len(data)) - if "true_h_max" not in data.colnames: # Avoid redundant warnings - self.log.warning("true_x_max not found, using 0") - - # Telescope position - # TableLoader with instrument=True usually adds pos_x, pos_y, pos_z in GroundFrame - if "pos_x" in data.colnames: - tel_x = data["pos_x"].quantity - tel_y = data["pos_y"].quantity - tel_z = data["pos_z"].quantity - else: - # Fallback to older names if necessary, but user suggested 'pos_x' - tel_x = data["tel_pos_x"].quantity - tel_y = data["tel_pos_y"].quantity - tel_z = data["tel_pos_z"].quantity - - # Tilted Frame Transformation - # correct core and telescope positions to be in the tilted frame - pointing_alt = data["telescope_pointing_altitude"].quantity - pointing_az = data["telescope_pointing_azimuth"].quantity - - # Create SkyCoord for pointing (one per event) - # TiltedGroundFrame needs a pointing direction. - - p_dir = AltAz(alt=pointing_alt, az=pointing_az) - tilted_frame = TiltedGroundFrame(pointing_direction=p_dir) - nominal_frame = NominalFrame(origin=p_dir) - - source_alt = data["true_alt"].quantity - source_az = data["true_az"].quantity - source_sky_coord = SkyCoord(alt=source_alt, az=source_az, frame=AltAz()) - source_nominal = source_sky_coord.transform_to(nominal_frame) - - source_x = source_nominal.fov_lon.to_value(u.radian) - source_y = source_nominal.fov_lat.to_value(u.radian) - - # Transform Core - # GroundFrame -> TiltedGroundFrame - # We assume core_z = 0 for the core position on ground? - # Usually true_core is on the ground plane (z=0 in GroundFrame). - core_ground = SkyCoord(x=core_x, y=core_y, z=0 * u.m, frame=GroundFrame()) - core_tilted = core_ground.transform_to(tilted_frame) - core_acc_x = core_tilted.x.to_value(u.m) - core_acc_y = core_tilted.y.to_value(u.m) - - # Transform Telescope Position - tel_ground = SkyCoord(x=tel_x, y=tel_y, z=tel_z, frame=GroundFrame()) - tel_tilted = tel_ground.transform_to(tilted_frame) - tel_acc_x = tel_tilted.x.to_value(u.m) - tel_acc_y = tel_tilted.y.to_value(u.m) - - # Calculate impact distance in Tilted Frame - dx = core_acc_x - tel_acc_x - dy = core_acc_y - tel_acc_y - phi = np.arctan2(dy, dx) - impact = np.sqrt(dx**2 + dy**2) - - # Memory Optimization: Use mask indices instead of tiling - # Get indices of valid pixels (rows=event_idx, cols=pixel_idx) - rows, cols = np.nonzero(mask_array) - - # Select features for valid pixels - energies_sel = energies[rows] - impact_sel = impact[rows] - xmax_sel = x_max[rows] - - # Select coordinates (Nominal Frame / FOV) - # fov_lon/lat are 1D arrays of pixel positions (Radians) - pix_lon_sel = fov_lon[cols] - pix_lat_sel = fov_lat[cols] - - # Select Source Position and Rotation Angle - # source_x/y and phi are per event (Radians) - source_x_sel = source_x[rows] - source_y_sel = source_y[rows] - phi_sel = phi[rows] - - # Rotation center on source - diff_lon = pix_lon_sel - source_x_sel - diff_lat = pix_lat_sel - source_y_sel - - # Rotate - cos_phi = np.cos(-phi_sel) - sin_phi = np.sin(-phi_sel) - - pix_lon_rot = diff_lon * cos_phi - diff_lat * sin_phi - pix_lat_rot = diff_lon * sin_phi + diff_lat * cos_phi - - # Select Images (Amplitudes) - images_sel = images[mask_array] - - # custom_symlog(np.array([images_sel/ 10]).ravel(), linear_threshold=2) - # tf.config.set_visible_devices([], "GPU") - # Stack features - X_chunk = np.stack( - [ - pix_lon_rot, - pix_lat_rot, - np.log10(energies_sel), - impact_sel / 100, - xmax_sel / 100, - custom_symlog( - np.array([images_sel / 10]), linear_threshold=2 - ).ravel(), - ], - axis=1, - ).astype(np.float32) - - pixel_features_list.append(X_chunk) - n_ev = np.sum(mask_array) - - n_loaded += n_ev if self.limit and n_loaded >= self.limit: break @@ -364,7 +469,6 @@ def _load_and_process_data(self, tel_type): return np.array([]), np.array([]) X_all = np.concatenate(pixel_features_list) - print(X_all.shape) # Dataset creation (Signal vs Background) n_total_pixels = len(X_all) n_samples = n_total_pixels // 2 @@ -389,27 +493,18 @@ def _load_and_process_data(self, tel_type): return X, y - def _build_model(self, input_shape): + def _build_model(self, input_shape, layer_number=10, layer_neurons=50): """Build Keras MLP.""" from tensorflow.keras.layers import Input - model = Sequential( - [ - Input(shape=input_shape), - Dense(64, activation="swish"), - Dense(64, activation="swish"), - Dense(64, activation="swish"), - Dense(64, activation="swish"), - Dense(64, activation="swish"), - Dense(64, activation="swish"), - Dense(64, activation="swish"), - Dense(64, activation="swish"), - Dense(64, activation="swish"), - Dense(64, activation="swish"), - Dense(64, activation="swish"), - Dense(1, activation="sigmoid"), - ] - ) + layers = [Input(shape=input_shape)] + + for _ in range(layer_number): + layers.append(Dense(layer_neurons, activation="swish")) + + layers.append(Dense(1, activation="sigmoid")) + + model = Sequential(layers) model.compile( optimizer=Adam(learning_rate=0.001), @@ -420,7 +515,7 @@ def _build_model(self, input_shape): def finish(self): """Cleanup.""" - self.loader.close() + pass def main(): From 5f766d63d15bc9870bcfaad9800ba08289603443 Mon Sep 17 00:00:00 2001 From: Dan Parsons Date: Thu, 5 Feb 2026 16:46:12 +0100 Subject: [PATCH 41/43] Add optional dependency check and fixed loading of subarray from hdf file --- src/ctapipe/tools/train_freepact.py | 31 +++++++++++++++++------------ 1 file changed, 18 insertions(+), 13 deletions(-) diff --git a/src/ctapipe/tools/train_freepact.py b/src/ctapipe/tools/train_freepact.py index 1f70db9484e..dd05eaf1179 100644 --- a/src/ctapipe/tools/train_freepact.py +++ b/src/ctapipe/tools/train_freepact.py @@ -4,14 +4,25 @@ import astropy.units as u import numpy as np -from tensorflow.keras.callbacks import EarlyStopping, ModelCheckpoint, ReduceLROnPlateau -from tensorflow.keras.layers import Dense -from tensorflow.keras.models import Sequential -from tensorflow.keras.optimizers import Adam + +try: + import tensorflow as tf + from tensorflow.keras.callbacks import ( + EarlyStopping, + ModelCheckpoint, + ReduceLROnPlateau, + ) + from tensorflow.keras.layers import Dense + from tensorflow.keras.models import Sequential + from tensorflow.keras.optimizers import Adam +except ImportError: + tf = None from ctapipe.core import Tool from ctapipe.core.traits import Float, Int, IntTelescopeParameter, List, Path +from ctapipe.exceptions import OptionalDependencyMissing from ctapipe.image import dilate, tailcuts_clean +from ctapipe.instrument import SubarrayDescription from ctapipe.io import TableLoader from ctapipe.utils.template_network_interpolator import custom_symlog @@ -29,6 +40,8 @@ class TrainFreePACT(Tool): name = "ctapipe-train-freepact" description = __doc__ + if tf is None: + raise OptionalDependencyMissing("tensorflow") output_path = Path( directory_ok=False, @@ -140,15 +153,7 @@ def setup(self): # Open the first file to get the subarray description # We assume all files have the same subarray try: - with TableLoader( - input_url=self.input_url[0], - parent=self, - dl1_images=True, - simulated=True, - true_parameters=True, - instrument=True, - ) as loader: - self.subarray = loader.subarray + self.subarray = SubarrayDescription.from_hdf(self.input_url[0]) except Exception as e: self.log.critical(f"Could not read subarray from first file: {e}") self.exit(1) From f1b81e79212ed0f594c84174daf6f41ca0894c4b Mon Sep 17 00:00:00 2001 From: Dan Parsons Date: Fri, 6 Feb 2026 14:40:29 +0100 Subject: [PATCH 42/43] Made tf import optional --- src/ctapipe/tools/train_freepact.py | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) diff --git a/src/ctapipe/tools/train_freepact.py b/src/ctapipe/tools/train_freepact.py index dd05eaf1179..f712ab500d5 100644 --- a/src/ctapipe/tools/train_freepact.py +++ b/src/ctapipe/tools/train_freepact.py @@ -40,8 +40,6 @@ class TrainFreePACT(Tool): name = "ctapipe-train-freepact" description = __doc__ - if tf is None: - raise OptionalDependencyMissing("tensorflow") output_path = Path( directory_ok=False, @@ -140,11 +138,8 @@ class TrainFreePACT(Tool): def setup(self): """Initialize components.""" # Check for tensorflow - try: - import tensorflow # noqa: F401 - except ImportError: - self.log.critical("tensorflow is required for this tool.") - self.exit(1) + if tf is None: + raise OptionalDependencyMissing("tensorflow") if not self.input_url: self.log.critical("Specifying input file(s) is required.") From b477ae556896dad336c8a4b4e8a7ddd92c2ee559 Mon Sep 17 00:00:00 2001 From: Dan Parsons Date: Tue, 24 Feb 2026 10:43:18 +0100 Subject: [PATCH 43/43] Some options added to FreePACT trainer, include trainer as command line executable --- pyproject.toml | 1 + src/ctapipe/tools/train_freepact.py | 79 +++++++++++++++++++++-------- 2 files changed, 59 insertions(+), 21 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index e46b27121de..120229183d0 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -120,6 +120,7 @@ ctapipe-calculate-pixel-statistics = "ctapipe.tools.calculate_pixel_stats:main" ctapipe-train-energy-regressor = "ctapipe.tools.train_energy_regressor:main" ctapipe-train-particle-classifier = "ctapipe.tools.train_particle_classifier:main" ctapipe-train-disp-reconstructor = "ctapipe.tools.train_disp_reconstructor:main" +ctapipe-train-freepact = "ctapipe.tools.train_freepact:main" ctapipe-apply-models = "ctapipe.tools.apply_models:main" ctapipe-store-astropy-cache = "ctapipe.tools.store_astropy_cache:main" diff --git a/src/ctapipe/tools/train_freepact.py b/src/ctapipe/tools/train_freepact.py index f712ab500d5..251d2308a5f 100644 --- a/src/ctapipe/tools/train_freepact.py +++ b/src/ctapipe/tools/train_freepact.py @@ -8,18 +8,20 @@ try: import tensorflow as tf from tensorflow.keras.callbacks import ( + CSVLogger, EarlyStopping, ModelCheckpoint, ReduceLROnPlateau, ) from tensorflow.keras.layers import Dense - from tensorflow.keras.models import Sequential + from tensorflow.keras.models import Sequential, load_model from tensorflow.keras.optimizers import Adam + except ImportError: tf = None from ctapipe.core import Tool -from ctapipe.core.traits import Float, Int, IntTelescopeParameter, List, Path +from ctapipe.core.traits import Bool, Float, Int, IntTelescopeParameter, List, Path from ctapipe.exceptions import OptionalDependencyMissing from ctapipe.image import dilate, tailcuts_clean from ctapipe.instrument import SubarrayDescription @@ -119,6 +121,17 @@ class TrainFreePACT(Tool): help="Additional rows to clean around the core.", ).tag(config=True) + max_threads = Int( + default_value=None, + allow_none=True, + help="Maximum number of CPU threads to use for training. Default is all available.", + ).tag(config=True) + + no_gpu = Bool( + default_value=False, + help="Do not use GPU for training, even if available.", + ).tag(config=True) + aliases = { ("i", "input"): "TrainFreePACT.input_url", ("o", "output"): "TrainFreePACT.output_path", @@ -131,6 +144,8 @@ class TrainFreePACT(Tool): "reweight_energy": "TrainFreePACT.reweight_energy", "layer_number": "TrainFreePACT.layer_number", "layer_neurons": "TrainFreePACT.layer_neurons", + "max_threads": "TrainFreePACT.max_threads", + "no_gpu": "TrainFreePACT.no_gpu", } classes = [TableLoader] @@ -141,6 +156,14 @@ def setup(self): if tf is None: raise OptionalDependencyMissing("tensorflow") + if self.max_threads is not None: + self.log.info("Setting maximum number of threads to %d", self.max_threads) + tf.config.threading.set_inter_op_parallelism_threads(self.max_threads) + + if self.no_gpu: + self.log.info("Disabling GPU usage") + tf.config.set_visible_devices([], "GPU") + if not self.input_url: self.log.critical("Specifying input file(s) is required.") self.exit(1) @@ -184,13 +207,17 @@ def start(self): input_shape=(X.shape[1],), layer_number=self.layer_number, layer_neurons=self.layer_neurons, + activation="swish", ) out_name = str(self.output_path).replace(".keras", "").replace(".h5", "") save_path = f"{out_name}_{tel_type}.keras" + log_path = save_path.replace(".keras", ".csv") + self.log.info("Training logs will be saved to %s", log_path) + callbacks = [ EarlyStopping( - monitor="val_loss", patience=50, verbose=1, start_from_epoch=50 + monitor="val_loss", patience=30, verbose=1, start_from_epoch=50 ), ReduceLROnPlateau( monitor="val_loss", @@ -205,7 +232,9 @@ def start(self): save_best_only=True, verbose=1, ), + CSVLogger(log_path, separator=",", append=False), ] + model.fit( X, y, @@ -215,6 +244,7 @@ def start(self): verbose=1, callbacks=callbacks, ) + model.layers[-1].activation = tf.keras.activations.linear # Save model # Construct a filename for this telescope type @@ -222,9 +252,13 @@ def start(self): # or save a dictionary of models? # Keras models are usually saved as directories or files. # Let's append the telescope type to the output name. + model.save(save_path) + model_temp = load_model(save_path) + + save_path = f"{out_name}_{tel_type}/" self.log.info("Saving model to %s", save_path) - model.save(save_path) + model_temp.export(save_path) def _load_and_process_data(self, tel_type): """Read chunked data and create pixel-wise dataset.""" @@ -263,6 +297,7 @@ def _load_and_process_data(self, tel_type): pixel_features_list = [] n_loaded = 0 + generator = np.random.default_rng() for input_file in self.input_url: self.log.info("Loading data from %s", input_file) @@ -305,7 +340,7 @@ def _load_and_process_data(self, tel_type): self.tailcut_threshold1, self.tailcut_threshold2, ) - for j in range(self.tailcut_additional_rows): + for _ in range(self.tailcut_additional_rows): mask = dilate(geometry, mask) mask_array[i] = mask @@ -314,7 +349,7 @@ def _load_and_process_data(self, tel_type): weight_factor = ( energies / (self.reweight_energy * u.TeV) ) ** self.reweight_index - random_value = np.random.rand(mask_array.shape[0]) + random_value = generator.random(mask_array.shape[0]) reweight_mask = random_value < weight_factor.value mask_array = np.logical_and( mask_array, reweight_mask[:, np.newaxis] @@ -468,39 +503,41 @@ def _load_and_process_data(self, tel_type): if not pixel_features_list: return np.array([]), np.array([]) - X_all = np.concatenate(pixel_features_list) + x_all = np.concatenate(pixel_features_list) # Dataset creation (Signal vs Background) - n_total_pixels = len(X_all) + n_total_pixels = len(x_all) n_samples = n_total_pixels // 2 - indices = np.random.choice(n_total_pixels, size=n_total_pixels, replace=False) + indices = generator.choice(n_total_pixels, size=n_total_pixels, replace=False) indices_sig = indices[:n_samples] indices_bg = indices[n_samples : 2 * n_samples] - X_sig = X_all[indices_sig] - y_sig = np.ones(len(X_sig)) + x_sig = x_all[indices_sig] + y_sig = np.ones(len(x_sig)) - X_bg = X_all[indices_bg].copy() - np.random.shuffle(X_bg[:, 5]) # Shuffle amplitude - y_bg = np.zeros(len(X_bg)) + x_bg = x_all[indices_bg].copy() + generator.shuffle(x_bg[:, 5]) # Shuffle amplitude + y_bg = np.zeros(len(x_bg)) - X = np.concatenate([X_sig, X_bg]) + x = np.concatenate([x_sig, x_bg]) y = np.concatenate([y_sig, y_bg]) - train_idx = np.random.permutation(len(X)) - X = X[train_idx] + train_idx = generator.permutation(len(x)) + x = x[train_idx] y = y[train_idx] - return X, y + return x, y - def _build_model(self, input_shape, layer_number=10, layer_neurons=50): + def _build_model( + self, input_shape, layer_number=10, layer_neurons=50, activation="swish" + ): """Build Keras MLP.""" from tensorflow.keras.layers import Input - layers = [Input(shape=input_shape)] + layers = [Input(shape=input_shape, name="input_1")] for _ in range(layer_number): - layers.append(Dense(layer_neurons, activation="swish")) + layers.append(Dense(layer_neurons, activation=activation)) layers.append(Dense(1, activation="sigmoid"))