Livetime scaling (#73)

allow changing livetime parameters by passing them to likelihood evaluations, data generation, fitting and confidence intervals

alea/models/blueice_extended_model.py

@@ -27,6 +27,8 @@ class BlueiceExtendedModel(StatisticalModel):
         is_data_set (bool): Whether data is set.
         _likelihood (LogLikelihoodSum): A blueice LogLikelihoodSum instance.
         likelihood_names (list): List of likelihood names.
+        livetime_parameter_names (list): List of the name of the livetime of each term,
+            None if not specified
         data_generators (list): List of data generators for each likelihood term.
 
     Args:
@@ -49,6 +51,9 @@ def __init__(self, parameter_definition: dict, likelihood_config: dict, **kwargs
         self._likelihood = self._build_ll_from_config(likelihood_config)
         self.likelihood_names = [t["name"] for t in likelihood_config["likelihood_terms"]]
         self.likelihood_names.append("ancillary_likelihood")
+        self.livetime_parameter_names = [t.get("livetime_parameter", None) for t in
+                                         likelihood_config["likelihood_terms"]]
+        self.livetime_parameter_names += [None]  # ancillary likelihood
         self.data_generators = self._build_data_generators()
 
     @classmethod
@@ -130,11 +135,14 @@ def get_expectation_values(self, **kwargs) -> dict:
         self_copy.data = self_copy.generate_data()
 
         # ancillary likelihood does not contribute
-        for ll_term, parameter_names in zip(
+        for ll_term, parameter_names, livetime_parameter in zip(  # noqa WPS352
                 self_copy._likelihood.likelihood_list[:-1],
-                self_copy._likelihood.likelihood_parameters):
+                self_copy._likelihood.likelihood_parameters,
+                self_copy.livetime_parameter_names):
             # WARNING: This silently drops parameters it can't handle!
             call_args = {k: i for k, i in generate_values.items() if k in parameter_names}
+            if livetime_parameter is not None:
+                call_args["livetime_days"] = generate_values[livetime_parameter]
 
             mus = ll_term(full_output=True, **call_args)[1]
             for n, mu in zip(ll_term.source_name_list, mus):
@@ -175,10 +183,13 @@ def _build_ll_from_config(self, likelihood_config: dict) -> "LogLikelihoodSum":
                     p.name for p in self.parameters if (
                         p.ptype == "shape") and (p.name not in source["parameters"])]
                 # no efficiency affects PDF:
-                parameters_to_ignore += [p.name for p in self.parameters if (p.ptype == "efficiency")]
+                parameters_to_ignore += [
+                    p.name for p in self.parameters if (
+                        p.ptype == "efficiency")]
                 parameters_to_ignore += source.get("extra_dont_hash_settings", [])
 
-                # ignore all shape parameters known to this model not named specifically in the source:
+                # ignore all shape parameters known to this model not named specifically
+                # in the source:
                 blueice_config["sources"][i]["extra_dont_hash_settings"] = parameters_to_ignore
 
             ll = likelihood_object(blueice_config)
@@ -239,7 +250,8 @@ def _build_data_generators(self) -> list:
             BlueiceDataGenerator(ll_term) for ll_term in self._likelihood.likelihood_list[:-1]]
 
     def _ll(self, **generate_values) -> float:
-        return self._likelihood(**generate_values)
+        livetime_days = [generate_values.get(ln, None) for ln in self.livetime_parameter_names]
+        return self._likelihood(livetime_days=livetime_days, **generate_values)
 
     def _generate_data(self, **generate_values) -> dict:
         """
@@ -270,10 +282,11 @@ def store_data(self, file_name, data_list, data_name_list=None, metadata=None):
             data_name_list = self.likelihood_names + ["generate_values"]
         super().store_data(file_name, data_list, data_name_list, metadata)
 
-    def _generate_science_data(self, **generate_values) -> dict:
+    def _generate_science_data(self,**generate_values)-> dict:
         """Generate the science data for all likelihood terms except the ancillary likelihood."""
-        science_data = [
-            gen.simulate(**generate_values) for gen in self.data_generators]
+        livetime_days = [generate_values.get(ln, None) for ln in self.livetime_parameter_names]
+        science_data = [gen.simulate(livetime_days=lt, **generate_values)
+                        for gen, lt in zip(self.data_generators, livetime_days)]
         return dict(zip(self.likelihood_names[:-1], science_data))
 
     def _generate_ancillary_measurements(self, **generate_values) -> dict:

alea/simulators.py

@@ -106,11 +106,12 @@ def simulate(
             The dtype follows self.dtype.
         """
         if filter_kwargs:
-            kwargs = {k: v for k, v in kwargs.items() if k in self.parameters}
+            kwargs = {k: v for k, v in kwargs.items() if k in self.parameters + ["livetime_days"]}
 
         unmatched_item = set(self.last_kwargs.items()) ^ set(kwargs.items())
         logging.debug("filtered kwargs in simulate: " + str(kwargs))
         logging.debug("unmatched_item in simulate: " + str(unmatched_item))
+        # check if the cached generator may be used:
         if ("FAKE_PARAMETER"
                 in self.last_kwargs.keys()) or (len(unmatched_item) != 0):
             ret = self.ll(full_output=True, **kwargs)  # result, mus, ps

alea/template_source.py

@@ -35,7 +35,7 @@ def build_histogram(self):
         histname = self.config['histname'].format(**format_dict)
 
         slice_args = self.config.get("slice_args", {})
-        if type(slice_args) is dict:
+        if isinstance(slice_args, dict):
             slice_args = [slice_args]
 
         h = template_to_multihist(templatename, histname)
@@ -57,7 +57,8 @@ def build_histogram(self):
                 slice_axis_limits = sa.get(
                     "slice_axis_limits", [bes[0], bes[-1]])
                 if sum_axis:
-                    logging.debug(f"Slice and sum over axis {slice_axis} from {slice_axis_limits[0]} to {slice_axis_limits[1]}")
+                    logging.debug(
+                        f"Slice and sum over axis {slice_axis} from {slice_axis_limits[0]} to {slice_axis_limits[1]}")
                     axis_names = h.axis_names_without(slice_axis)
                     h = h.slicesum(
                         axis=slice_axis,
@@ -66,7 +67,8 @@ def build_histogram(self):
                     h.axis_names = axis_names
                 else:
                     logging.debug(f"Normalization before slicing: {h.n}.")
-                    logging.debug(f"Slice over axis {slice_axis} from {slice_axis_limits[0]} to {slice_axis_limits[1]}")
+                    logging.debug(
+                        f"Slice over axis {slice_axis} from {slice_axis_limits[0]} to {slice_axis_limits[1]}")
                     h = h.slice(axis=slice_axis,
                                 start=slice_axis_limits[0],
                                 stop=slice_axis_limits[1])
@@ -126,7 +128,8 @@ def build_histogram(self):
         self._n_events_histogram = h.similar_blank_histogram()
 
         h *= self.config.get('histogram_scale_factor', 1)
-        logging.debug(f"Multiplying histogram with histogram_scale_factor {self.config.get('histogram_scale_factor', 1)}. Histogram is now normalised to {h.n}.")
+        logging.debug(
+            f"Multiplying histogram with histogram_scale_factor {self.config.get('histogram_scale_factor', 1)}. Histogram is now normalised to {h.n}.")
 
         # Convert h to density...
         if self.config.get('in_events_per_bin'):
@@ -147,7 +150,7 @@ def build_histogram(self):
         if np.min(self._pdf_histogram.histogram) < 0:
             raise AssertionError(
                 f"There are bins for source {templatename} with negative entries."
-                    )
+            )
 
     def simulate(self, n_events):
         dtype = []
@@ -189,9 +192,9 @@ def build_histogram(self):
         templatenames = self.config['templatenames']
 
         slice_args = self.config.get("slice_args", {})
-        if type(slice_args) is dict:
+        if isinstance(slice_args, dict):
             slice_args = [slice_args]
-        if type(slice_args[0]) is dict:
+        if isinstance(slice_args[0], dict):
             slice_argss = []
             for n in histnames:
                 slice_argss.append(slice_args)
@@ -387,11 +390,11 @@ def build_histogram(self):
         histname = self.config['histname'].format(**format_dict)
 
         spectrum = self.config["spectrum"]
-        if type(spectrum) is str:
+        if isinstance(spectrum, str):
             spectrum = self._get_json_spectrum(spectrum.format(**format_dict))
 
         slice_args = self.config.get("slice_args", {})
-        if type(slice_args) is dict:
+        if isinstance(slice_args, dict):
             slice_args = [slice_args]
 
         h = template_to_multihist(templatename, histname)