diff --git a/configs/cascade_7param_MCPE.yaml b/configs/cascade_7param_MCPE.yaml new file mode 100644 index 00000000..4d2b5e8f --- /dev/null +++ b/configs/cascade_7param_MCPE.yaml @@ -0,0 +1,426 @@ +--- +############## +# Config for single cascade events +# +# Detailed information on the parameters are given in the SetupManager +# class located in egenerator.settings.setup_manager.py. +############## + +# Provide a unique name for the model +unique_name: 'cascade_7param_MCPE' + + +#------------------ +# Training settings +#------------------ +training_settings: { + 'optimizer_name': 'Adam', + 'optimizer_settings': { + # 'amsgrad': True, + 'learning_rate': { + 'full_class_string': 'egenerator.utils.learning_rate.MultiLearningRateScheduler', + 'settings':{ + 'boundaries': [10000, 500000], + 'scheduler_settings': [ + { + 'full_class_string': 'tensorflow.keras.optimizers.schedules.PolynomialDecay', + 'settings': { + 'initial_learning_rate': 0.01, + 'end_learning_rate': 0.001, + 'decay_steps': 10000, + }, + }, + { + 'full_class_string': 'tensorflow.keras.optimizers.schedules.PolynomialDecay', + 'settings': { + 'initial_learning_rate': 0.001, + 'end_learning_rate': 0.001, + 'decay_steps': 490000, + }, + }, + { + 'full_class_string': 'tensorflow.keras.optimizers.schedules.PolynomialDecay', + 'settings': { + 'initial_learning_rate': 0.001, + 'end_learning_rate': 0.000001, + 'decay_steps': 500000, + 'power': 2, + }, + }, + ] + }, + }, + }, + 'l1_regularization': 0., + 'l2_regularization': 0., + 'clip_gradients_value': , + 'remove_nan_gradients': False, + + 'validation_frequency': 100, + 'evaluation_frequency': 100, + 'save_frequency': 500, + + # Number of training iterations to train the model for + 'num_training_iterations': 1000000, + + # Additional keywords to the loss module used for training + 'additional_loss_module_kwargs': { + 'normalize_by_total_charge': False, + }, +} + +#----------------.------- +# Reconstruction settings +#------------------------ +reconstruction_settings: { + 'reco_output_file': '../data/reconstruction/{unique_name}/reconstruction_{unique_name}.hdf5', + + # Define which modules to run + 'calculate_covariance_matrix': False, + 'calculate_goodness_of_fit': False, + 'estimate_angular_uncertainty': False, + 'run_mcmc': False, + 'make_1d_llh_scans': False, + + # define which tensor to use as parameter_tensor (default: x_parameters) + 'parameter_tensor_name': 'x_parameters', + # define seed (to seed from MC truth choose: 'x_parameters') + 'seed': 'MonopodFit4_PartialExclusion', + # these are only relevant for the SeedLoaderMiscModule. This specifies + # from which column names in the data table the parameters will be loaded + 'seed_parameter_names': ['x', 'y', 'z', 'zenith', 'azimuth', + 'energy', 'time', + ], + 'seed_missing_value': 1., + 'seed_missing_value_dict': { + }, + 'seed_float_precision': 'float32', + + # choose the otpimizer iterface: + # 'scipy' or 'tfp' (tensorflow_probability) + 'reco_optimizer_interface': 'scipy', + + 'scipy_optimizer_settings': { + # 'method': 'L-BFGS-B', + 'method': 'BFGS', + 'options': { + # 'ftol': !!float 1e-7, + } + }, + + 'tf_optimizer_settings': { + 'method': 'bfgs_minimize', + 'x_tolerance': 0.001, + }, + + # Perform minimization in transformed and normalized parameter space + # if True. This is usually desired as it may facilitate minimization. + 'minimize_in_trafo_space': True, + + # Specify which parameters to fit. + # If True, the parameter is minimized doing reconstruction, otherwise it + # is held constant. + # Set default value which will apply to all parameters, except if stated + # otherwise in the 'minimize_parameter_dict' + 'minimize_parameter_default_value': True, + # Settings defined here overwrite the default value + # Entries must have the form: {parameter_name: value} + 'minimize_parameter_dict': { + + }, +} + +#----------------------- +# Model Manager settings +#----------------------- + +# Settings for model manager class +model_manager_settings: { + + # The loss module class to use + 'model_manager_class': 'egenerator.manager.source.SourceManager', + + # restore model if True, otherwise start from scratch + 'restore_model' : True, + + # These settings are used to configure the model manager and may not change + config: { + + # Path to where the manager will be saved to and loaded from + 'manager_dir': '../data/training/{unique_name}/manager', + + # Define which model parameters to use for the I3Particle + # [x, y, z, zenith, azimuth, energy, time] + 'I3ParticleMapping': { + 'x': 'cascade_x', + 'y': 'cascade_y', + 'z': 'cascade_z', + 'zenith': 'cascade_zenith', + 'azimuth': 'cascade_azimuth', + 'energy': 'cascade_energy', + 'time': 'cascade_time', + }, + }, +} + +#----------------------- +# Data Iterator settings +#----------------------- + +# These settings describe the data iterators +data_iterator_settings: { + + # The training data iterator + 'training': { + 'batch_size': 32, + 'num_splits': , + 'file_capacity': 10, + 'batch_capacity': 500, + 'dataset_capacity': 500, + 'num_jobs': 10, + 'num_add_files': 5, + 'num_repetitions': 3, + 'pick_random_files_forever': False, + 'input_data': [ + '/data/user/fschluter/eg-data/electrons/unbinned/set_3/Gen1k_10TeV_EMinus_FTP_V3_*.hdf5', + '/data/user/fschluter/eg-data/electrons/unbinned/set_4/Gen1k_10TeV_EMinus_FTP_V3_*.hdf5', + '/data/user/fschluter/eg-data/electrons/unbinned/set_5/Gen1k_10TeV_EMinus_FTP_V3_*.hdf5', + '/data/user/fschluter/eg-data/electrons/unbinned/set_6/Gen1k_10TeV_EMinus_FTP_V3_*.hdf5', + '/data/user/fschluter/eg-data/electrons/unbinned/set_7/Gen1k_10TeV_EMinus_FTP_V3_*.hdf5', + '/data/user/fschluter/eg-data/electrons/unbinned/set_8/Gen1k_10TeV_EMinus_FTP_V3_*.hdf5', + '/data/user/fschluter/eg-data/electrons/unbinned/set_9/Gen1k_10TeV_EMinus_FTP_V3_*.hdf5', + '/data/user/fschluter/eg-data/electrons/unbinned/set_10/Gen1k_10TeV_EMinus_FTP_V3_*.hdf5', + '/data/user/fschluter/eg-data/electrons/unbinned/set_11/Gen1k_10TeV_EMinus_FTP_V3_*.hdf5', + '/data/user/fschluter/eg-data/electrons/unbinned/set_12/Gen1k_10TeV_EMinus_FTP_V3_*.hdf5', + '/data/user/fschluter/eg-data/electrons/unbinned/set_13/Gen1k_10TeV_EMinus_FTP_V3_*.hdf5', + '/data/user/fschluter/eg-data/electrons/unbinned/set_14/Gen1k_10TeV_EMinus_FTP_V3_*.hdf5', + '/data/user/fschluter/eg-data/electrons/unbinned/set_15/Gen1k_10TeV_EMinus_FTP_V3_*.hdf5', + ], + }, + + # The trafo data iterator + 'trafo': { + 'batch_size': 32, + 'num_splits': , + 'file_capacity': 1, + 'batch_capacity': 2, + 'num_jobs': 12, + 'num_add_files': 1, + 'num_repetitions': 1, + 'pick_random_files_forever': False, + 'input_data': [ + '/data/user/fschluter/eg-data/electrons/unbinned/set_3/Gen1k_10TeV_EMinus_FTP_V3_*.hdf5', + '/data/user/fschluter/eg-data/electrons/unbinned/set_4/Gen1k_10TeV_EMinus_FTP_V3_*.hdf5', + '/data/user/fschluter/eg-data/electrons/unbinned/set_5/Gen1k_10TeV_EMinus_FTP_V3_*.hdf5', + '/data/user/fschluter/eg-data/electrons/unbinned/set_6/Gen1k_10TeV_EMinus_FTP_V3_*.hdf5', + '/data/user/fschluter/eg-data/electrons/unbinned/set_7/Gen1k_10TeV_EMinus_FTP_V3_*.hdf5', + '/data/user/fschluter/eg-data/electrons/unbinned/set_8/Gen1k_10TeV_EMinus_FTP_V3_*.hdf5', + '/data/user/fschluter/eg-data/electrons/unbinned/set_9/Gen1k_10TeV_EMinus_FTP_V3_*.hdf5', + '/data/user/fschluter/eg-data/electrons/unbinned/set_10/Gen1k_10TeV_EMinus_FTP_V3_*.hdf5', + '/data/user/fschluter/eg-data/electrons/unbinned/set_11/Gen1k_10TeV_EMinus_FTP_V3_*.hdf5', + '/data/user/fschluter/eg-data/electrons/unbinned/set_12/Gen1k_10TeV_EMinus_FTP_V3_*.hdf5', + '/data/user/fschluter/eg-data/electrons/unbinned/set_13/Gen1k_10TeV_EMinus_FTP_V3_*.hdf5', + '/data/user/fschluter/eg-data/electrons/unbinned/set_14/Gen1k_10TeV_EMinus_FTP_V3_*.hdf5', + '/data/user/fschluter/eg-data/electrons/unbinned/set_15/Gen1k_10TeV_EMinus_FTP_V3_*.hdf5', + ], + }, + + # The validation data iterator + 'validation': { + 'batch_size': 32, + 'num_splits': , + 'file_capacity': 10, + 'batch_capacity': 10, + 'dataset_capacity': 10, + 'num_jobs': 10, + 'num_add_files': 5, + 'num_repetitions': 1, + 'pick_random_files_forever': False, + 'input_data': [ + '/data/user/fschluter/eg-data/electrons/unbinned/set_0/Gen1k_10TeV_EMinus_FTP_V3_*.hdf5', + '/data/user/fschluter/eg-data/electrons/unbinned/set_1/Gen1k_10TeV_EMinus_FTP_V3_*.hdf5', + '/data/user/fschluter/eg-data/electrons/unbinned/set_2/Gen1k_10TeV_EMinus_FTP_V3_*.hdf5', + ], + }, + + # The test data iterator + 'test': { + 'batch_size': 1, + 'num_splits': , + 'file_capacity': 1, + 'batch_capacity': 10, + 'num_jobs': 1, + 'num_add_files': 0, + 'num_repetitions': 1, + 'sample_randomly': False, + 'pick_random_files_forever': False, + 'input_data': [ + '/data/user/fschluter/eg-data/electrons/unbinned/set_16/Gen1k_10TeV_EMinus_FTP_V3_*.hdf5', + '/data/user/fschluter/eg-data/electrons/unbinned/set_17/Gen1k_10TeV_EMinus_FTP_V3_*.hdf5', + '/data/user/fschluter/eg-data/electrons/unbinned/set_18/Gen1k_10TeV_EMinus_FTP_V3_*.hdf5', + ], + }, +} + +#--------------------- +# Loss module settings +#--------------------- + +# Settings for the loss module class +loss_module_settings: { + # The loss module class to use + 'loss_module': 'egenerator.loss.default.DefaultLossModule', + + config: { + # the float precision to use + 'float_precision': 'float32', + + 'add_normalization_term': True, + + # choose the loss function to use + 'loss_function_name': 'unbinned_extended_pulse_llh', # this implies a Poisson likelihood for the charge + }, +} + +#--------------------------- +# Evaluation module settings +#--------------------------- + +# Settings for the evaluation module class +evaluation_module_settings: { + + # The loss module class to use + 'evaluation_module':, + config: { + }, +} + +#----------------------------- +# Data Transformation settings +#----------------------------- +data_trafo_settings: { + 'float_precision': 'float64', + 'norm_constant': !!float 1e-6, + 'num_batches': 2000, + 'model_dir': '../data/trafo/{unique_name}', +} + +#---------------------- +# Data Handler settings +#---------------------- + +# Settings for the data handler class +data_handler_settings: { + + # The data handler class to use + 'data_handler': 'modular.ModuleDataHandler', + + # ------------------- + # DataModule Settings + # ------------------- + # which data module to use + 'data_module': 'pulse_data.PulseDataModule', + + # settings for the data module + 'data_settings':{ + 'pulse_key': 'I3MCPESeriesMap', + 'pulse_is_mcpe': True, + 'dom_exclusions_key':, + 'time_exclusions_key':, + 'float_precision': 'float32', + 'add_charge_quantiles': False, + 'discard_pulses_from_excluded_doms': False, + }, + + # -------------------- + # LabelModule Settings + # -------------------- + # which label module to use + 'label_module': 'cascades.CascadeGeneratorLabelModule', + + # settings for the label module + 'label_settings':{ + 'shift_cascade_vertex': True, + # logarithm on labels: + # (x, y, z, zenith, azimuth, energy, time)? + # 'parameter_names': ['x', 'y', 'z', 'zenith', 'azimuth', 'energy', 'time'], + 'trafo_log': [False, False, False, False, False, True, False], + 'label_key': 'LabelsDeepLearning', + 'float_precision': 'float32', + }, + + # --------------------- + # WeightModule Settings + # --------------------- + # which weight module to use + 'weight_module': 'dummy.DummyWeightModule', + + # settings for the weight module + 'weight_settings':{}, + + # ------------------- + # MiscModule Settings + # ------------------- + # which misc module to use + 'misc_module': 'dummy.DummyMiscModule', + + # settings for the misc module + 'misc_settings':{}, + + # ------------------- + # FilterModule Settings + # ------------------- + # which filter module to use + 'filter_module': 'general_filter.GeneralFilterModule', + + # settings for the filter module + 'filter_settings':{ + 'constraints': [ + # ['LabelsDeepLearning', 'cascade_z', '<', -250.], + ], + }, +} + +#--------------- +# Model settings +#--------------- + +# Settings for the neural network model class +model_settings: { + + # The source class to use + 'model_class': 'egenerator.model.source.cascade.default.DefaultCascadeModel', + + config: { + 'keep_prob':, + 'add_opening_angle': True, + 'add_dom_coordinates': False, + 'num_local_vars': 0, + 'scale_charge': True, + 'scale_charge_by_relative_dom_efficiency': True, + 'scale_charge_by_global_dom_efficiency': False, + 'prevent_mixture_component_swapping': False, + 'estimate_charge_distribution': 'negative_binomial', + 'num_latent_models': 10, + + # This is a list of labels in addition to + # (x, y, z, zenith, azimuth, energy, time) and snowstorm parameters + 'additional_label_names' : [ + ], + + # First convolutions + 'filter_size_list' : [[1, 1], [1, 1], [1, 1], [1, 1], [1, 1], [1, 1]], + 'num_filters_list' : [25, 500, 500, 500, 42], # The last dimension should be 4 * num_asym_gauss + num_charge_pdf; num_asym_gauss = 10, num_charge_pdf = 1 (poisson) or 2 (binominal) + 'method_list' : ['locally_connected', + 'convolution', 'convolution', 'convolution', + 'convolution', + ], + 'use_dropout_list' : False, + 'activation_list' : ['elu', 'elu', 'elu', 'elu', ''], + 'use_batch_norm_list' : False, + 'use_residual_list' : True, + }, +} + +#---------------------- +... diff --git a/egenerator/data/handler/modular.py b/egenerator/data/handler/modular.py index b5b514fa..9e91aab6 100644 --- a/egenerator/data/handler/modular.py +++ b/egenerator/data/handler/modular.py @@ -343,8 +343,10 @@ def _get_data(self, file_or_frame, method, *args, **kwargs): assert isinstance(file_or_frame, str), "Expected file path string" num_data, data = self.data_module.get_data_from_hdf( - file_or_frame, *args, **kwargs - ) + file_or_frame, *args, + label_key=self.label_module.configuration.config["label_key"], + **kwargs) + num_labels, labels = self.label_module.get_data_from_hdf( file_or_frame, *args, **kwargs ) diff --git a/egenerator/data/modules/data/pulse_data.py b/egenerator/data/modules/data/pulse_data.py index 0081f789..739c9f98 100644 --- a/egenerator/data/modules/data/pulse_data.py +++ b/egenerator/data/modules/data/pulse_data.py @@ -45,6 +45,7 @@ def _configure( float_precision, add_charge_quantiles, discard_pulses_from_excluded_doms, + pulse_is_mcpe=False, ): """Configure Module Class This is an abstract method and must be implemented by derived class. @@ -72,6 +73,9 @@ def _configure( discard_pulses_from_excluded_doms : bool, optional If True, pulses on excluded DOMs are discarded. The pulses are discarded after the charge at the DOM is collected. + pulse_is_mcpe : bool, optional + If True, train on MCPE pulses instead of RecoPulses. This setting + changes how the charge is read out of the hdf5 files. (Default: False) Returns ------- @@ -214,6 +218,7 @@ def _configure( ), mutable_settings=dict( pulse_key=pulse_key, + pulse_is_mcpe=pulse_is_mcpe, dom_exclusions_key=dom_exclusions_key, time_exclusions_key=time_exclusions_key, discard_pulses_from_excluded_doms=( @@ -247,12 +252,13 @@ def get_data_from_hdf(self, file, *args, **kwargs): if not self.is_configured: raise ValueError("Module not configured yet!") - # open file - f = pd.HDFStore(file, "r") + charge_str = "npe" if self.configuration.config["pulse_is_mcpe"] else "charge" + # open file + f = pd.HDFStore(file, 'r') try: pulses = f[self.configuration.config["pulse_key"]] - _labels = f["LabelsDeepLearning"] + _labels = f[kwargs["label_key"]] if self.data["dom_exclusions_exist"]: try: dom_exclusions = f[ @@ -293,11 +299,13 @@ def get_data_from_hdf(self, file, *args, **kwargs): # create Dictionary with event IDs size = len(_labels["Event"]) + + if not size: + raise ValueError("Label length is 0.") + event_dict = {} - for idx, row in _labels.iterrows(): - event_dict[ - (row.iloc[0], row.iloc[1], row.iloc[2], row.iloc[3]) - ] = idx + for row in _labels.itertuples(): + event_dict[(row[1:5])] = row[0] # create empty array for DOM charges x_dom_charge = np.zeros( @@ -351,21 +359,22 @@ def get_data_from_hdf(self, file, *args, **kwargs): "skipping pulse: {} {}".format(string, dom) ) continue + index = event_dict[(row[1:5])] # accumulate charge in DOMs - x_dom_charge[index, string - 1, dom - 1, 0] += row.charge + x_dom_charge[index, string - 1, dom - 1, 0] += getattr(row, charge_str) # gather pulses if add_charge_quantiles: # (charge, time, quantile) cum_charge = float(x_dom_charge[index, string - 1, dom - 1, 0]) - x_pulses[pulse_index] = [row.charge, row.time, cum_charge] + x_pulses[pulse_index] = [getattr(row, charge_str), row.time, cum_charge] else: # (charge, time) - x_pulses[pulse_index] = [row.charge, row.time] + x_pulses[pulse_index] = [getattr(row, charge_str), row.time] # gather pulse ids (batch index, string, dom) x_pulses_ids[pulse_index] = [index, string - 1, dom - 1] @@ -408,7 +417,7 @@ def get_data_from_hdf(self, file, *args, **kwargs): continue index = event_dict[(row[1:5])] - # t_start (pulse time): row[10], t_end (pulse width): row[11] + # t_start (pulse time): row.time, t_end (pulse width): row[11] # (t_start, t_end) x_time_exclusions[tw_index] = [row.time, row.width] @@ -577,7 +586,6 @@ def get_data_from_frame(self, frame, *args, **kwargs): for pulse in pulse_list: index = 0 - # pulse charge: row[12], time: row[10] # accumulate charge in DOMs x_dom_charge[index, string - 1, dom - 1, 0] += pulse.charge diff --git a/egenerator/data/modules/labels/cascades.py b/egenerator/data/modules/labels/cascades.py index b3554470..cce1b3a4 100644 --- a/egenerator/data/modules/labels/cascades.py +++ b/egenerator/data/modules/labels/cascades.py @@ -34,6 +34,9 @@ def _configure( trafo_log, float_precision, label_key="LabelsDeepLearning", + parameter_names=["cascade_x", "cascade_y", "cascade_z", + "cascade_zenith", "cascade_azimuth", + "cascade_energy", "cascade_t"], ): """Configure Module Class This is an abstract method and must be implemented by derived class. @@ -51,10 +54,12 @@ def _configure( If a single bool is given, this applies to all labels. Otherwise a list of bools corresponds to the labels in the order: x, y, z, zenith, azimuth, energy, time - label_key : str, optional - The name of the key under which the labels are saved. float_precision : str The float precision as a str. + label_key : str, optional + The name of the key under which the labels are saved. + parameter_names : list of str, optional + Name of the parameters (e.g, the columns in the hdf5 Dataset `label_key`) Returns ------- @@ -129,6 +134,7 @@ def _configure( trafo_log=trafo_log, float_precision=float_precision, label_key=label_key, + parameter_names=parameter_names, ), ) return configuration, data, {} @@ -163,16 +169,8 @@ def get_data_from_hdf(self, file, *args, **kwargs): cascade_parameters = [] try: _labels = f[self.configuration.config["label_key"]] - for label in [ - "cascade_x", - "cascade_y", - "cascade_z", - "cascade_zenith", - "cascade_azimuth", - "cascade_energy", - "cascade_t", - ]: - cascade_parameters.append(_labels[label]) + for par in self.configuration.config["parameter_names"]: + cascade_parameters.append(_labels[par]) except Exception as e: self._logger.warning(e) @@ -223,16 +221,8 @@ def get_data_from_frame(self, frame, *args, **kwargs): cascade_parameters = [] try: _labels = frame[self.configuration.config["label_key"]] - for label in [ - "cascade_x", - "cascade_y", - "cascade_z", - "cascade_zenith", - "cascade_azimuth", - "cascade_energy", - "cascade_t", - ]: - cascade_parameters.append(np.atleast_1d(_labels[label])) + for par in self.configuration.config["parameter_names"]: + cascade_parameters.append(np.atleast_1d(_labels[par])) except Exception as e: self._logger.warning(e) diff --git a/egenerator/ic3/__init__.py b/egenerator/ic3/__init__.py index 2228e2ef..b56d37b5 100644 --- a/egenerator/ic3/__init__.py +++ b/egenerator/ic3/__init__.py @@ -1,7 +1,9 @@ from egenerator.ic3.reconstruction import EventGeneratorReconstruction from egenerator.ic3.simulation import EventGeneratorSimulation +from egenerator.ic3.evaluate_mcpe import CalculateLikelihood __all__ = [ "EventGeneratorReconstruction", "EventGeneratorSimulation", + "CalculateLikelihood", ] diff --git a/egenerator/ic3/evaluate_mcpe.py b/egenerator/ic3/evaluate_mcpe.py new file mode 100644 index 00000000..99dbae24 --- /dev/null +++ b/egenerator/ic3/evaluate_mcpe.py @@ -0,0 +1,621 @@ +import os +import numpy as np +import tensorflow as tf +import timeit +import math + + +import matplotlib.pyplot as plt + +from icecube import dataclasses, icetray, photonics_service, simclasses +from icecube.icetray.i3logging import log_info, log_warn + +from egenerator.utils.configurator import ManagerConfigurator +from egenerator.utils import basis_functions +from egenerator.ic3.simulation import EventGeneratorSimulation + +def calculate_loss(pred, true): + assert true % 1 == 0, "Not a int value" + return np.exp(-1 * pred) * pred ** true / math.factorial(true) + + +class CalculateLikelihood(EventGeneratorSimulation): + + """Class to evaluate predictions of MCPE with Event-Generator model. + + """ + + def __init__(self, context): + """Class to evaluate predictions of MCPE with Event-Generator model. + + Parameters + ---------- + context : TYPE + Description + """ + icetray.I3ConditionalModule.__init__(self, context) + self.AddOutBox('OutBox') + + # Required settings + self.AddParameter('model_name', + 'A model name that defines which model to ' + 'apply. The full path is given by `model_base_dir` +' + ' model_name if `model_base_dir` is provided. ' + 'Otherwise `model_name` must define the full ' + 'model path.') + + # Optional settings + self.AddParameter('default_values', + 'Not all parameters of the source hypothesis can ' + 'be extracted from the I3MCTree. For these ' + 'parameters, default values may be defined. The ' + '`default_values` must be a dictionary of the ' + 'format `parameter_name`: `value`. `value` may ' + 'either be a double or a string. If a string is ' + 'provided, it is assumed that an I3Double exists in ' + 'the frame under the key as provided by `value`', + { + 'Absorption': 1., + 'AnisotropyScale': 1., + 'DOMEfficiency': 1., + 'HoleIceForward_Unified_00': 0., + 'HoleIceForward_Unified_01': 0., + 'Scattering': 1., + }) + self.AddParameter('model_base_dir', + 'The base directory in which the model is located.' + ' The full path is given by `model_base_dir` + ' + ' model_name if `model_base_dir` is provided. ' + 'Otherwise `model_name` must define the full ' + 'model path.', + '/data/user/mhuennefeld/exported_models/egenerator') + self.AddParameter('output_key', + 'The output base key to which pulses will be saved', + 'EventGeneratorPrediction') + self.AddParameter('mc_tree_name', + 'The name of the propagated I3MCTree for which the ' + 'light yield and measured pulses will be simulated', + 'I3MCTree') + self.AddParameter('max_simulation_distance', + 'Particles further away from the center of the ' + 'detector than this distance are not predicted', + 2250) + self.AddParameter('min_simulation_energy', + 'Particles below this energy threshold are not ' + 'predicted.', + 100) + self.AddParameter('num_threads', + 'Number of threads to use for tensorflow. This will ' + 'be passed on to tensorflow where it is used to set ' + 'the "intra_op_parallelism_threads" and ' + '"inter_op_parallelism_threads" settings. If a ' + 'value of zero (default) is provided, the system ' + 'uses an appropriate number. Note: when running ' + 'this as a job on a cluster you might want to limit ' + '"num_threads" to the amount of allocated CPUs.', + 0) + self.AddParameter('random_service', + 'The random service or seed to use. If this is an ' + 'integer, a numpy random state will be created with ' + 'the seed set to `random_service`', + 42) + self.AddParameter('SimulateElectronDaughterParticles', + 'If true, look for daughter particles of an electron and ' + 'simulate their light yield rather than from the ' + 'mother electron. This is necessary to correctly ' + 'simulate the cascade longitudinal extension of ' + 'high energy cascades.', + False) + self.AddParameter('AddDarkDOMs', '', True) + self.AddParameter('ValidateTimePDF', '', True) + + self.AddParameter('CascadeService', + '', + None) + + + def Configure(self): + """Configures Module and loads model from file. + """ + self.model_name = self.GetParameter('model_name') + self.default_values = self.GetParameter('default_values') + self.model_base_dir = self.GetParameter('model_base_dir') + self.output_key = self.GetParameter('output_key') + self.mc_tree_name = self.GetParameter('mc_tree_name') + self.max_simulation_distance = \ + self.GetParameter('max_simulation_distance') + self.min_simulation_energy = self.GetParameter('min_simulation_energy') + self.num_threads = self.GetParameter('num_threads') + self.random_service = self.GetParameter('random_service') + self.simulate_electron_daughters = self.GetParameter('SimulateElectronDaughterParticles') + self.add_dark_doms = self.GetParameter("AddDarkDOMs") + self.validate_time_pdf = self.GetParameter("ValidateTimePDF") + if isinstance(self.random_service, int): + self.random_service = np.random.RandomState(self.random_service) + + if self.model_base_dir is not None: + self.model_dir = os.path.join(self.model_base_dir, self.model_name) + else: + self.model_dir = self.model_name + + # -------------------------------------------------- + # Build and configure SourceManager and extrac Model + # -------------------------------------------------- + self.manager_configurator = ManagerConfigurator( + manager_dirs=[self.model_dir], + num_threads=self.num_threads, + ) + self.manager = self.manager_configurator.manager + + for model in self.manager.models: + num_vars, num_total_vars = model.num_variables + msg = '\nNumber of Model Variables:\n' + msg += '\tFree: {}\n' + msg += '\tTotal: {}' + log_info(msg.format(num_vars, num_total_vars)) + + if len(self.manager.models) > 1: + raise NotImplementedError( + 'Currently does not support model ensemble.') + + self.model = self.manager.models[0] + + # get parameter names that have to be set + self.param_names = sorted([n for n in self.model.parameter_names + if n not in self.default_values]) + + # make sure that the only parameters that need to be set are provided + included_parameters = [ + 'azimuth', 'energy', 'time', 'x', 'y', 'z', 'zenith'] + + # search if there is a common prefix that we can use + prefix_list = [] + for name in included_parameters: + if name in self.param_names: + prefix_list.append('') + else: + found_match = False + for param in self.param_names: + if param[-len(name):] == name: + # only allow prefixes ending on '_' + if param[-len(name)-1] == '_': + prefix_list.append(param[:-len(name)]) + found_match = True + if not found_match: + msg = ( + 'Did not find a parameter name match for "{}". Model ' + 'Parameter names are: {}' + ).format(name, self.param_names) + raise ValueError(msg) + + prefix_list = np.unique(prefix_list) + if len(prefix_list) != 1: + msg = 'Could not find common parameter prefix. Found: {}'.format( + prefix_list) + raise ValueError(msg) + + self._prefix = prefix_list[0] + + # double check that the prefix now works + for name in self.param_names: + if name[len(self._prefix):] not in included_parameters: + raise KeyError('Unknown parameter name:', name) + + # Create concrete tensorflow function to obtain DOM expectations + self.param_dtype = getattr( + tf, self.manager.data_trafo.data['tensors']['x_parameters'].dtype) + self.get_model_tensors = self.manager.get_model_tensors_function() + + # --------------------------------------------------- + # Define which particles are simulated by which Model + # --------------------------------------------------- + + # define allowed parent particles + # These are particles which we can safely simulate by only looking at + # their daughters + self.allowed_parent_particles = [ + dataclasses.I3Particle.NuE, + dataclasses.I3Particle.NuMu, + dataclasses.I3Particle.NuTau, + dataclasses.I3Particle.NuEBar, + dataclasses.I3Particle.NuMuBar, + dataclasses.I3Particle.NuTauBar, + dataclasses.I3Particle.Hadrons, + ] + if self.simulate_electron_daughters: + self.allowed_parent_particles += [ + dataclasses.I3Particle.EMinus, + dataclasses.I3Particle.EPlus, + ] + + # Define type of particles that can be simulated as EM cascades + self.em_cascades = [ + dataclasses.I3Particle.Pi0, + dataclasses.I3Particle.Gamma, + dataclasses.I3Particle.PairProd, + dataclasses.I3Particle.Brems, + dataclasses.I3Particle.DeltaE, + dataclasses.I3Particle.EMinus, + dataclasses.I3Particle.EPlus, + ] + + # Define type of particles that can be simulated as tracks + self.tracks = [ + dataclasses.I3Particle.MuMinus, + dataclasses.I3Particle.MuPlus, + ] + + # define particles that do not deposit light, in other words + # if we end up with a final state particle of this type it is ok + # not to simulate the light yield for these particles. + self.dark_particles = [ + dataclasses.I3Particle.NuE, + dataclasses.I3Particle.NuMu, + dataclasses.I3Particle.NuTau, + dataclasses.I3Particle.NuEBar, + dataclasses.I3Particle.NuMuBar, + dataclasses.I3Particle.NuTauBar, + ] + + self.cascade_service = self.GetParameter('CascadeService') + + if self.cascade_service is not None: + self.photonics = Photonics(self.cascade_service) + + self._counter = 0 + + def DAQ(self, frame): + """Apply Event-Generator model to physics frames. + + Parameters + ---------- + frame : I3Frame + The current Q-Frame. + """ + + self._counter += 1 + + if not self._counter % 100: + print(self._counter) + + # start timer + t_0 = timeit.default_timer() + + # get sources + cascades, tracks = self.get_light_sources(frame[self.mc_tree_name]) + + if len(tracks) > 0: + raise NotImplementedError('Tracks not yet supported') + + # convert cascade to source hypotheses + cascade_sources = self.convert_cascades_to_tensor(cascades, frame) + + # timer after source collection + t_1 = timeit.default_timer() + + # get result tensors from model + result_tensors = self.get_model_tensors(cascade_sources) + + # timer after NN evaluation + t_2 = timeit.default_timer() + + # sample pdf + if self.validate_time_pdf: + self.sample_from_pdf(frame, result_tensors, validate_photonics=self.cascade_service is not None) + if self._counter % 10 == 0: + self.plot_pdfs(frame, result_tensors, validate_photonics=self.cascade_service is not None) + else: + self.store_total_charge(frame, result_tensors, validate_photonics=self.cascade_service is not None) + + # timer after Sampling + t_3 = timeit.default_timer() + + log_info('Simulation took: {:3.3f}ms'.format((t_3 - t_0) * 1000)) + log_info('\t Gathering Sources: {:3.3f}ms'.format((t_1 - t_0) * 1000)) + log_info('\t Evaluating NN model: {:3.3f}ms'.format((t_2 - t_1) * 1000)) + log_info('\t Sampling Pulses: {:3.3f}ms'.format((t_3 - t_2) * 1000)) + + # push frame to next modules + self.PushFrame(frame) + + def store_total_charge(self, frame, result_tensors, validate_photonics=False): + """ This function stores the predicted charge by EG, MC, Photonics """ + total_charge = result_tensors['dom_charges'].numpy() + mcpe_series_map = frame["I3MCPESeriesMap"] + + # We have to use here I3MCPulseSeriesMap here because for a simple I3MapKeyDouble is no hdf5 converter implemented ... + + event_generator_total_dom_charge = dataclasses.I3MapKeyVectorDouble() + mc_total_dom_charge = dataclasses.I3MapKeyVectorDouble() + + if validate_photonics: + geo = frame["I3Geometry"] + particle = frame["I3MCTree"].get_head() + photonics_total_dom_charge = dataclasses.I3MapKeyVectorDouble() + + for string in range(86): + for om in range(60): + omkey = icetray.OMKey(string + 1, om + 1) + event_generator_total_dom_charge[omkey] = [float(total_charge[0, string, om, 0])] + + if validate_photonics: + position = geo.omgeo[omkey].position + self.photonics.cascade_pxs.SelectModuleCoordinates(*position) + charge, _, t0 = self.photonics.cascade_pxs.SelectSource(photonics_service.PhotonicsSource(particle)) + photonics_total_dom_charge[omkey] = [float(charge)] + + if omkey in mcpe_series_map: + mc_total_dom_charge[omkey] = [np.sum([pulse.npe for pulse in mcpe_series_map[omkey]], dtype="float")] + else: + mc_total_dom_charge[omkey] = [0] + + frame.Put("MCTotalChargePerDOM", mc_total_dom_charge) + frame.Put("EventGeneratorTotalChargePerDOM", event_generator_total_dom_charge) + + if validate_photonics: + frame.Put("PhotonicsTotalChargePerDOM", photonics_total_dom_charge) + + + def sample_from_pdf(self, frame, result_tensors, validate_photonics=False): + total_charge = result_tensors['dom_charges'].numpy() + eps = 1e-7 + mcpe_series_map = frame["I3MCPESeriesMap"] + + loss_event_generator = dataclasses.I3MapKeyDouble() + loss_event_generator_total = 0 + + event_generator_total_dom_charge = dataclasses.I3MapKeyDouble() + mc_total_dom_charge = dataclasses.I3MapKeyDouble() + + if validate_photonics: + loss_photonics = dataclasses.I3MapKeyDouble() + loss_photonics_total = 0 + geo = frame["I3Geometry"] + particle = frame["I3MCTree"].get_head() + photonics_total_dom_charge = dataclasses.I3MapKeyDouble() + + for string in range(86): + for om in range(60): + omkey = icetray.OMKey(string + 1, om + 1) + + event_generator_total_dom_charge[omkey] = float(total_charge[0, string, om, 0]) + + if omkey in mcpe_series_map: + pulse_series = mcpe_series_map[omkey] + + dom_charges_true = np.array([pulse.npe for pulse in pulse_series]) + mc_total_dom_charge[omkey] = float(np.sum(dom_charges_true)) + pulse_log_pdf_values = np.log([self.model.pdf_per_dom(pulse.time, result_tensors=result_tensors, + string=string, dom=om) + eps for pulse in pulse_series]) + + dom_charges_pred = total_charge[0, string, om, 0] + + time_log_likelihood = -dom_charges_true * pulse_log_pdf_values + llh_poisson = dom_charges_pred - np.sum(dom_charges_true) * np.log(dom_charges_pred + eps) + + loss = np.sum(time_log_likelihood) + np.sum(llh_poisson) + + loss_event_generator_total += loss + loss_event_generator[omkey] = loss + + if validate_photonics: + position = geo.omgeo[omkey].position + # for pulse in pulse_series: + table_pdf, dom_charges_pred = self.photonics.get_pdf( + particle, position, [pulse.time for pulse in pulse_series], + quantiles=False) + photonics_total_dom_charge[omkey] = float(dom_charges_pred) + + time_log_likelihood = -dom_charges_true * np.log(np.array(table_pdf) + eps) + llh_poisson = dom_charges_pred - np.sum(dom_charges_true) * np.log(dom_charges_pred + eps) + + loss = np.sum(time_log_likelihood) + np.sum(llh_poisson) + + loss_photonics_total += loss + loss_photonics[omkey] = loss + + elif self.add_dark_doms: + expected_charge = total_charge[0, string, om, 0] + mc_total_dom_charge[omkey] = 0 + loss = expected_charge - 0 * np.log(expected_charge + 1e-7) + loss_event_generator_total += loss + loss_event_generator[omkey] = loss + + if validate_photonics: + position = geo.omgeo[omkey].position + expected_charge, _ = self.photonics.get_PE(particle, position) + photonics_total_dom_charge[omkey] = float(expected_charge) + + loss = expected_charge - 0 * np.log(expected_charge + 1e-7) + loss_photonics_total += loss + loss_photonics[omkey] = loss + + else: + pass + + frame.Put("MCTotalChargePerDOM", mc_total_dom_charge) + frame.Put("EventGeneratorTotalChargePerDOM", event_generator_total_dom_charge) + frame.Put("EventGeneratorLoss", loss_event_generator) + frame.Put("EventGeneratorLossTotal", dataclasses.I3Double(loss_event_generator_total)) + + if validate_photonics: + frame.Put("PhotonicsTotalChargePerDOM", photonics_total_dom_charge) + frame.Put("PhotonicsLoss", loss_photonics) + frame.Put("PhotonicsLossTotal", dataclasses.I3Double(loss_photonics_total)) + + + def plot_pdfs(self, frame, result_tensors, validate_photonics=False): + + total_charge = result_tensors['dom_charges'].numpy() + mcpe_series_map = frame["I3MCPESeriesMap"] + + if validate_photonics: + geo = frame["I3Geometry"] + particle = frame["I3MCTree"].get_head() + + hom = None + max_charge = 0 + + for omkey, pulses in mcpe_series_map: + tot_charge = np.sum(p.npe for p in pulses) + + if tot_charge > max_charge: + max_charge = tot_charge + hom = omkey + + if hom is None: + return + + string = hom.string - 1 + om = hom.om - 1 + + if validate_photonics: + geo = frame["I3Geometry"] + particle = frame["I3MCTree"].get_head() + position = geo.omgeo[hom].position + self.photonics.cascade_pxs.SelectModuleCoordinates(*position) + _, _, t0 = self.photonics.cascade_pxs.SelectSource(photonics_service.PhotonicsSource(particle)) + else: + t0 = 0 + + times = np.arange(t0 - 50, t0 + 300, 5) + dt = np.diff(times)[0] + ctime = times[:-1] + dt / 2 + + fig, ax = plt.subplots() + charges_mc = [p.npe for p in mcpe_series_map[hom]] + ax.hist([p.time for p in mcpe_series_map[hom]], times, + weights=charges_mc, alpha=0.5, color="k", label=f"MC {np.sum(charges_mc)}") + + eg_pdf_values = np.squeeze(self.model.get_probability_quantiles_per_dom(times, result_tensors=result_tensors, + string=string, dom=om)) + tot_charge = total_charge[0, string, om, 0] + ax.plot(ctime, eg_pdf_values * tot_charge, label=f"EventGenerator {tot_charge:.2f}", lw=1) + + if validate_photonics: + + table_pdf, dom_charges_pred = self.photonics.get_pdf(particle, position, times) + ax.plot(ctime, np.array(table_pdf) * dom_charges_pred, label=f"Photonics {dom_charges_pred:.2f}", lw=1, ls=":") + + + ax.set_xlabel("time / ns") + ax.set_ylabel("npe") + ax.legend() + ax.grid() + + header = frame["I3EventHeader"] + + plt.tight_layout() + plt.savefig(f"{self.model_name}_{header.run_id}_{header.event_id}_pdf.png") + + +class Photonics(object): + """ + Little wrapper around photonics_service + """ + + def __init__(self, cascadePhotonicsService=None): + + if cascadePhotonicsService is None: + + if os.environ.get('I3_DATA') is None or os.environ.get('I3_DATA') == "": + raise ValueError("The environment variable \"I3_DATA\" " + "is not set. Are you in an IceTray environment?") + + table_base = os.path.expandvars('$I3_DATA/photon-tables/splines/ems_spice1_z20_a10.%s.fits') + cascadePhotonicsService = photonics_service.I3PhotoSplineService( + table_base % 'abs', table_base % 'prob', timingSigma=0) + + self.cascade_pxs = cascadePhotonicsService + + + def _get_PE(self, source, position): + """ Returns the number of photons and geometrical distance for a source & reciever pair """ + + if not isinstance(source, dataclasses.I3Particle): + source = convert_params_to_particle(source) + + if source.energy == 0: + return 0, 0 + + self.cascade_pxs.SelectModuleCoordinates(*position) + pes, dist, _ = self.cascade_pxs.SelectSource(photonics_service.PhotonicsSource(source)) + if pes < 0: + pes = 0 + + return pes, dist + + + def get_PE(self, sources, position): + """ Wrapper around self._get_PE to allow providing a list of sources """ + + if isinstance(sources, dataclasses.I3Particle): + return self._get_PE(sources, position) + + elif isinstance(sources, list) and isinstance(sources[0], float): + return self._get_PE(sources, position) + + elif isinstance(sources, list): + output = np.array([self._get_PE(source, position) for source in sources]) + pe_sum = np.sum(output[:, 0]) + mask = output[:, 0] > 0 + if not np.any(mask): + return 0, 0 + return pe_sum, np.amin(output[mask, 1]) + else: + raise TypeError(f"The type of sources \"{type(sources)}\" is not supported.") + + + def _get_pdf(self, source, position, times, quantiles=True): + + if not isinstance(source, dataclasses.I3Particle): + source = convert_params_to_particle(source) + + if source.energy == 0: + return np.zeros_like(times), 0 + + self.cascade_pxs.SelectModuleCoordinates(*position) + pes, _, t0 = self.cascade_pxs.SelectSource(photonics_service.PhotonicsSource(source)) + if pes <= 0: + return np.zeros_like(times), 0 + + if quantiles: + pdf = self.cascade_pxs.GetProbabilityQuantiles(times, t0, 0) # Reading charge distribution (normalized to 1) + + # dt = np.diff(times) + # ctime = times[:-1] + dt / 2 + # pdf2 = [self.cascade_pxs.GetProbabilityDensity(float(t) - t0) for t in ctime] + # pdf2 = pdf2 * dt + + # for t, p1, p2 in zip(ctime, pdf, pdf2): + # print(f"{t:.1f}, {p1 / p2:.2f}") + # print() + + return pdf, pes + else: + if isinstance(times, list): + pdf = [self.cascade_pxs.GetProbabilityDensity(t - t0) for t in times] + else: + pdf = self.cascade_pxs.GetProbabilityDensity(times - t0) + return pdf, pes + + + def get_pdf(self, sources, position, t, quantiles=True): + if isinstance(sources, dataclasses.I3Particle): + return self._get_pdf(sources, position, t, quantiles) + + elif isinstance(sources, list) and isinstance(sources[0], float): + return self._get_pdf(sources, position, t, quantiles) + + else: + pdf_tot = np.zeros_like(t) + pdfs = [] + nps_tot = 0 + for source in sources: + pdf, nps = self._self._get_pdf(source, position, t, quantiles) + pdfs.append(pdf) + pdf_tot += pdf * nps + nps_tot += nps + + pdfs.append(pdf_tot / nps_tot) + return pdfs \ No newline at end of file diff --git a/egenerator/ic3/reconstruction.py b/egenerator/ic3/reconstruction.py index 4984b32d..7afb2b2b 100644 --- a/egenerator/ic3/reconstruction.py +++ b/egenerator/ic3/reconstruction.py @@ -389,12 +389,7 @@ def Configure(self): self.manager = self.manager_configurator.manager self.loss_module = self.manager_configurator.loss_module - if "I3ParticleMapping" in self.manager.configuration.config["config"]: - self.i3_mapping = self.manager.configuration.config["config"][ - "I3ParticleMapping" - ] - else: - self.i3_mapping = None + self.i3_mapping = self.manager.configuration.config["config"].get("I3ParticleMapping", None) for model in self.manager.models: num_vars, num_total_vars = model.num_variables diff --git a/egenerator/ic3/visualization.py b/egenerator/ic3/visualization.py index 704e4065..11a6ae21 100644 --- a/egenerator/ic3/visualization.py +++ b/egenerator/ic3/visualization.py @@ -128,12 +128,7 @@ def Configure(self): self.manager = self.manager_configurator.manager self.loss_module = self.manager_configurator.loss_module - if "I3ParticleMapping" in self.manager.configuration.config["config"]: - self.i3_mapping = self.manager.configuration.config["config"][ - "I3ParticleMapping" - ] - else: - self.i3_mapping = None + self.i3_mapping = self.manager.configuration.config["config"].get("I3ParticleMapping", None) for model in self.manager.models: num_vars, num_total_vars = model.num_variables diff --git a/egenerator/manager/base.py b/egenerator/manager/base.py index 18e1afbf..e7577a66 100644 --- a/egenerator/manager/base.py +++ b/egenerator/manager/base.py @@ -641,10 +641,7 @@ def perform_training_step( gradients = tape.gradient(combined_loss, variables) # remove nans in gradients and replace these with zeros - if "remove_nan_gradients" in opt_config: - remove_nan_gradients = opt_config["remove_nan_gradients"] - else: - remove_nan_gradients = False + remove_nan_gradients = opt_config.get("remove_nan_gradients", False) if remove_nan_gradients: gradients = [ tf.where(tf.is_nan(grad), tf.zeros_like(grad), grad) @@ -652,10 +649,7 @@ def perform_training_step( if grad is not None ] - if "clip_gradients_value" in opt_config: - clip_gradients_value = opt_config["clip_gradients_value"] - else: - clip_gradients_value = None + clip_gradients_value = opt_config.get("clip_gradients_value", None) if clip_gradients_value is not None: capped_gradients, _ = tf.clip_by_global_norm( gradients, clip_gradients_value @@ -675,6 +669,7 @@ def perform_training_step( ], ) asserts.append(assert_finite) + with tf.control_dependencies(asserts): self.optimizer.apply_gradients(zip(capped_gradients, variables)) @@ -778,10 +773,8 @@ def train( scheduler_class = misc.load_class(lr_cfg["full_class_string"]) scheduler = scheduler_class(**lr_cfg["settings"]) optimizer_settings["learning_rate"] = scheduler - - optimizer = getattr(tf.optimizers, opt_config["optimizer_name"])( - **optimizer_settings - ) + + optimizer = getattr(tf.optimizers, opt_config["optimizer_name"])(**optimizer_settings) self._untracked_data["optimizer"] = optimizer # save new training step to model @@ -866,7 +859,6 @@ def train( # increment step counter for model in self.models: model.step.assign_add(1) - # get new batch of training data training_data_batch = next(train_dataset) @@ -876,7 +868,7 @@ def train( # -------------------------- # evaluate on validation set # -------------------------- - if step % opt_config["validation_frequency"] == 0: + if step % opt_config["validation_frequency"] == 0 and step: # not validate on the fist iteration new_validation_time = timeit.default_timer() time_diff = new_validation_time - validation_time validation_time = new_validation_time diff --git a/egenerator/model/base.py b/egenerator/model/base.py index 6c97da05..e0dbdb25 100644 --- a/egenerator/model/base.py +++ b/egenerator/model/base.py @@ -251,7 +251,7 @@ def save_weights( in the output directory. If it does not exist yet, a new one will be created. Otherwise, its values will be updated The file contains meta data on the checkpoints and keeps track - of the most recents files. The structure and content of meta data: + of the most recents files. The structure and content of meta data: latest_checkpoint: int The number of the latest checkpoint. diff --git a/egenerator/model/source/base.py b/egenerator/model/source/base.py index d6234cc8..def821c9 100644 --- a/egenerator/model/source/base.py +++ b/egenerator/model/source/base.py @@ -544,6 +544,132 @@ def cdf( return cdf_values + def cdf_per_dom(self, x, result_tensors, string, dom, **kwargs): + """Compute CDF values at x for given result_tensors + + This is a numpy, i.e. not tensorflow, method to compute the CDF based + on a provided `result_tensors`. This can be used to investigate + the generated PDFs. + + Note: this function only works for sources that use asymmetric + Gaussians to parameterize the PDF. The latent values of the AG + must be included in the `result_tensors`. + + Parameters + ---------- + x : array_like + The times in ns at which to evaluate the result tensors. + Shape: () or [n_points] + result_tensors : dict of tf.tensor + The dictionary of output tensors as obtained from `get_tensors`. + string : int + String id of the DOM + dom : int + Module id of the DOM + **kwargs + Keyword arguments. + + Returns + ------- + array_like + The CDF values at times x for the given event hypothesis and + exclusions that were used to compute `result_tensors`. + Shape: [n_events, n_points] + + Raises + ------ + NotImplementedError + If assymetric Gaussian latent variables are not present in + `result_tensors` dictionary. + """ + + x_orig = np.atleast_1d(x) + assert len(x_orig.shape) == 1, x_orig.shape + n_points = len(x_orig) + + # shape: [1, 1, n_points] + x = np.reshape(x_orig, (1, 1, -1)) + + if result_tensors["time_offsets"] is not None: + # shape: [n_events] + t_offsets = result_tensors["time_offsets"].numpy() + # shape: [n_events, 1, 1] + t_offsets = np.reshape(t_offsets, [-1, 1, 1]) + # shape: [n_events, 1, n_points] + x = x - t_offsets + else: + t_offsets = 0. + + # internally we are working with different time units + x = x / self.time_unit_in_ns + + # Check if the asymmetric Gaussian latent variables exist + for latent_name in ["mu", "scale", "sigma", "r"]: + if "latent_var_" + latent_name not in result_tensors: + msg = "PDF evaluation is not supported for this model: {}" + raise NotImplementedError(msg.format( + self._untracked_data["name"])) + + # extract values + # shape: [n_events, n_components, 1] + mu = result_tensors["latent_var_mu"].numpy()[:, string, dom, :, np.newaxis] + scale = result_tensors["latent_var_scale"].numpy()[:, string, dom, :, np.newaxis] + sigma = result_tensors["latent_var_sigma"].numpy()[:, string, dom, :, np.newaxis] + r = result_tensors["latent_var_r"].numpy()[:, string, dom, :, np.newaxis] + + # shape: [n_events, n_components, n_points] + mixture_cdf = basis_functions.asymmetric_gauss_cdf( + x=x, mu=mu, sigma=sigma, r=r + ) + + # shape: [n_events, n_points] + cdf_values = np.sum(mixture_cdf * scale, axis=1) + + return cdf_values + + + def get_probability_quantiles_per_dom(self, x, result_tensors, string, dom, **kwargs): + """ + This is a numpy, i.e. not tensorflow, method to compute the CDF based + on a provided `result_tensors`. This can be used to investigate + the generated PDFs. + + Note: this function only works for sources that use asymmetric + Gaussians to parameterize the PDF. The latent values of the AG + must be included in the `result_tensors`. + + Parameters + ---------- + x : array_like + The times in ns at which to evaluate the result tensors. + Shape: () or [n_points] + result_tensors : dict of tf.tensor + The dictionary of output tensors as obtained from `get_tensors`. + string : int + String id of the DOM + dom : int + Module id of the DOM + **kwargs + Keyword arguments. + + Returns + ------- + array_like + Shape: [n_events, n_points - 1] + + Raises + ------ + NotImplementedError + If assymetric Gaussian latent variables are not present in + `result_tensors` dictionary. + """ + + cdf = self.cdf_per_dom(x, result_tensors, string, dom, **kwargs) + quantiles = cdf[..., 1:] - cdf[..., :-1] + + return quantiles + + def pdf( self, x, @@ -674,6 +800,95 @@ def pdf( return pdf_values + def pdf_per_dom(self, x, result_tensors, string, dom, **kwargs): + """Compute PDF values at x for given result_tensors for a specific dom + + This is a numpy, i.e. not tensorflow, method to compute the PDF based + on a provided `result_tensors`. This can be used to investigate + the generated PDFs. + + Note: this function only works for sources that use asymmetric + Gaussians to parameterize the PDF. The latent values of the AG + must be included in the `result_tensors`. + + Note: the PDF does not set values inside excluded time windows to zero, + but it does adjust the normalization. It is assumed that pulses will + already be masked before evaluated by Event-Generator. Therefore, an + extra check for exclusions is not performed due to performance issues. + + Parameters + ---------- + x : array_like + The times in ns at which to evaluate the result tensors. + Shape: () or [n_points] + result_tensors : dict of tf.tensor + The dictionary of output tensors as obtained from `get_tensors`. + string : int + String id of the DOM + dom : int + Module id of the DOM + **kwargs + Keyword arguments. + + Returns + ------- + array_like + The PDF values at times x for the given event hypothesis and + exclusions that were used to compute `result_tensors`. + Shape: [n_events, n_points] + + Raises + ------ + NotImplementedError + If assymetric Gaussian latent variables are not present in + `result_tensors` dictionary. + """ + + x_orig = np.atleast_1d(x) + assert len(x_orig.shape) == 1, x_orig.shape + + # shape: [1, 1, n_points] + x = np.reshape(x_orig, (1, 1, -1)) + + if result_tensors["time_offsets"] is not None: + # shape: [n_events] + t_offsets = result_tensors["time_offsets"].numpy() + + # shape: [n_events, 1, 1] + t_offsets = np.reshape(t_offsets, [-1, 1, 1]) + # shape: [n_events, 1, n_points] + x = x - t_offsets + + # internally we are working with different time units + x = x / self.time_unit_in_ns + + # Check if the asymmetric Gaussian latent variables exist + for latent_name in ["mu", "scale", "sigma", "r"]: + if "latent_var_" + latent_name not in result_tensors: + msg = "PDF evaluation is not supported for this model: {}" + raise NotImplementedError(msg.format( + self._untracked_data["name"])) + + # extract values + # shape: [n_events, n_components, 1] + mu = result_tensors["latent_var_mu"].numpy()[:, string, dom, :, np.newaxis] + scale = result_tensors["latent_var_scale"].numpy()[:, string, dom, :, np.newaxis] + sigma = result_tensors["latent_var_sigma"].numpy()[:, string, dom, :, np.newaxis] + r = result_tensors["latent_var_r"].numpy()[:, string, dom, :, np.newaxis] + + # shape: [n_events, n_components, n_points] + mixture_pdf = basis_functions.asymmetric_gauss( + x=x, mu=mu, sigma=sigma, r=r + ) + + # shape: [n_events, n_points] + pdf_values = np.sum(mixture_pdf * scale, axis=1) + + # invert back to PDF in ns + pdf_values = pdf_values / self.time_unit_in_ns + + return np.squeeze(pdf_values) + def _apply_pdf_time_window_exclusions( self, times, pdf_values, tw_exclusions, tw_exclusions_ids ): @@ -731,7 +946,7 @@ def _apply_pdf_time_window_exclusions( # Description # """ # parent_node = tf.constant(np.ones(shape=[1, self.num_parameters])) - # output = self._untracked_data['module'](parent_node) + # output = self._untracked_data["module"](parent_node) # return _find_top_nodes(output) # def _find_top_nodes(self, tensor): diff --git a/egenerator/model/source/cascade/default.py b/egenerator/model/source/cascade/default.py index 2cff3c09..b15543ce 100644 --- a/egenerator/model/source/cascade/default.py +++ b/egenerator/model/source/cascade/default.py @@ -378,12 +378,11 @@ def get_tensors( # check if we have the right amount of filters in the latent dimension n_models = config["num_latent_models"] - if n_models * 4 + n_charge != config["num_filters_list"][-1]: - raise ValueError( - "{!r} != {!r}".format( - n_models * 4 + n_charge, config["num_filters_list"][-1] - ) - ) + exp_num_filts = n_models * 4 + n_charge + if exp_num_filts != config["num_filters_list"][-1]: + raise ValueError("{!r} (expected) != {!r} (configured)".format( + exp_num_filts, config["num_filters_list"][-1])) + if n_models <= 1: raise ValueError("{!r} !> 1".format(n_models)) diff --git a/egenerator/utils/learning_rate.py b/egenerator/utils/learning_rate.py index 5e2062bc..9ab84ced 100644 --- a/egenerator/utils/learning_rate.py +++ b/egenerator/utils/learning_rate.py @@ -58,7 +58,6 @@ def __init__(self, boundaries, scheduler_settings, name=None): self.schedulers = schedulers self.name = name - @tf.function def __call__(self, step): step = tf.convert_to_tensor(step) @@ -86,7 +85,10 @@ def __call__(self, step): pred_fn_pairs.append( ( - step > low and step <= high, + tf.math.logical_and( + tf.math.greater(step, low), + tf.math.less_equal(step, high), + ), lambda: scheduler(step - low), ) )