Shape parameters

.gitignore

@@ -12,9 +12,11 @@ __pycache__
 *.pdf
 *.ipynb
 !notebooks/*.ipynb
+!alea/examples/*.ipynb
 *.h5
 *.hdf5
-!templates/*.h5
-!templates/*.hdf5
+!alea/examples/templates/*.h5
+!alea/examples/templates/*.hdf5
 docs/build
 build
+debug.py

alea/__init__.py

@@ -1,15 +1,7 @@
 __version__ = '0.0.0'
 
-from . import model
 from .model import *
-
 from . import models
-
-from . import utils
 from .utils import *
-
-from . import parameters
 from .parameters import *
-
-from . import simulators
 from .simulators import *

alea/_likelihoods/ll_nt_from_config.py

@@ -9,7 +9,7 @@
 import numpy as np
 import scipy.stats as sps
 from alea._plotting import pdf_plotter
-from alea.simulators import simulate_interpolated
+from alea.simulators import BlueiceDataGenerator
 from blueice.inference import bestfit_scipy, one_parameter_interval
 from blueice.likelihood import LogLikelihoodSum
 from inference_interface import (dict_to_structured_array,
@@ -272,13 +272,13 @@ def set_toy_reading_mode(self, toydata_mode, toydata_file):
         self.toydata_mode = toydata_mode
         self.toydata_file = toydata_file
         if toydata_mode == "none":
-            self.rgs = [simulate_interpolated(ll, binned=b) for ll, b in zip(self.lls, self.binned)]
+            self.rgs = [BlueiceDataGenerator(ll, binned=b) for ll, b in zip(self.lls, self.binned)]
             if hasattr(self, "datasets_array"):
                 del self.datasets_array
             if hasattr(self, "toydata_index"):
                 del self.toydata_index
         elif toydata_mode == "write":
-            self.rgs = [simulate_interpolated(ll, binned=b) for ll, b in zip(self.lls, self.binned)]
+            self.rgs = [BlueiceDataGenerator(ll, binned=b) for ll, b in zip(self.lls, self.binned)]
             self.datasets_array = []
             if hasattr(self, "toydata_index"):
                 del self.toydata_index

alea/examples/__init__.py

@@ -0,0 +1 @@
+from .gaussian_model import *

alea/model_configs/unbinned_wimp_statistical_model.yaml → alea/examples/model_configs/unbinned_wimp_statistical_model.yaml

@@ -45,24 +45,26 @@ parameter_definition:
     fit_guess: 1.0
     description: Parameter to account for the uncertain signal expectation given a certain cross-section
 
-  # er_band_shift:
-  #   nominal_value: 0
-  #   ptype: shape
-  #   uncertainty: scipy.stats.uniform(loc=-1, scale=2)
-  #   relative_uncertainty: false
-  #   fittable: true
-  #   blueice_anchors:
-  #     - -1
-  #     - 0
-  #     - 1
-  #   fit_limits:
-  #     - -1
-  #     - 1
-  #   description: ER band shape parameter (shifts the ER band up and down)
+  er_band_shift:
+    nominal_value: 0
+    ptype: shape
+    uncertainty: null  # scipy.stats.uniform(loc=-1, scale=2)
+    # relative_uncertainty: false
+    fittable: true
+    blueice_anchors:
+      - -2
+      - -1
+      - 0
+      - 1
+      - 2
+    fit_limits:
+      - -2
+      - 2
+    description: ER band shape parameter (shifts the ER band up and down)
 
 likelihood_config:
   likelihood_weights: [1, 1, 1]
-  template_folder: alea/templates
+  template_folder: alea/examples/templates
   likelihood_terms:
     # SR0
     - name: sr0
@@ -79,8 +81,10 @@ likelihood_config:
         histname: er_template  # TODO: implement a default histname based on the source name
         parameters:
           - er_rate_multiplier
-          # - er_band_shift
-        template_filename: er_template.h5
+          - er_band_shift
+        named_parameters:
+          - er_band_shift
+        template_filename: er_template_{er_band_shift}.h5
         histogram_scale_factor: 1
 
       - name: wimp
@@ -108,8 +112,10 @@ likelihood_config:
         histname: er_template
         parameters:
           - er_rate_multiplier
-          # - er_band_shift
-        template_filename: er_template.h5
+          - er_band_shift
+        named_parameters:
+          - er_band_shift
+        template_filename: er_template_{er_band_shift}.h5
         histogram_scale_factor: 2
 
       - name: wimp

alea/model.py

@@ -6,7 +6,6 @@
 from scipy.stats import chi2
 from scipy.optimize import brentq
 from iminuit import Minuit
-from iminuit.util import make_func_code
 from inference_interface import toydata_to_file
 
 from alea.parameters import Parameters
@@ -230,28 +229,27 @@ class MinuitWrap:
             the dictionary passed to the Minuit call.
             """
 
-            def __init__(self, f, s_args):
+            def __init__(self, f, parameters: Parameters):
                 self.func = f
-                self.s_args = s_args
-                self.func_code = make_func_code(s_args)
+                self.s_args = parameters.names
+                self._parameters = {p.name: p.fit_limits for p in parameters}
 
             def __call__(self, *args):
                 return self.func(args)
 
         # Make the Minuit object
-        cost.errordef = Minuit.LIKELIHOOD
         m = Minuit(
-            MinuitWrap(cost, s_args=self.parameters.names),
+            MinuitWrap(cost, parameters=self.parameters),
             **defaults)
+        m.errordef = Minuit.LIKELIHOOD
         fixed_params = [] if fixed_parameters is None else fixed_parameters
         fixed_params += self.parameters.not_fittable
         for par in fixed_params:
             m.fixed[par] = True
-        for n, l in self.parameters.fit_limits.items():
-            m.limits[n] = l
 
         # Call migrad to do the actual minimization
         m.migrad()
+        self.minuit_object = m
         if verbose:
             print(m)
         # alert! This gives the _maximum_ likelihood

alea/models/__init__.py

@@ -1,5 +1 @@
-from . import gaussian_model
-from .gaussian_model import *
-
-from . import blueice_extended_model
 from .blueice_extended_model import *

alea/models/blueice_extended_model.py

@@ -5,7 +5,7 @@
 import numpy as np
 import scipy.stats as stats
 from blueice.likelihood import LogAncillaryLikelihood, LogLikelihoodSum
-from inference_interface import dict_to_structured_array
+from inference_interface import dict_to_structured_array, structured_array_to_dict
 
 from alea.model import StatisticalModel
 from alea.simulators import BlueiceDataGenerator
@@ -56,21 +56,24 @@ def from_config(cls, config_file_path: str) -> "BlueiceExtendedModel":
         return cls(**config)
 
     @property
-    def data(self) -> list:
+    def data(self) -> dict:
         """Returns the data of the statistical model."""
         return super().data
 
     @data.setter
-    def data(self, data: list):
+    def data(self, data: dict):
         """
         Overrides default setter. Will also set the data of the blueice ll.
         Data-sets are expected to be in the form of a list of one
         or more structured arrays-- representing the data-sets of one or more likelihood terms.
         """
         # iterate through all likelihood terms and set the science data in the blueice ll
         # last entry in data are the generate_values
-        for d, ll_term in zip(data[:-1], self._likelihood.likelihood_list):
-            ll_term.set_data(d)
+        for i, (dataset_name, d) in enumerate(data.items()):
+            if dataset_name != "generate_values":
+                ll_term = self._likelihood.likelihood_list[i]
+                assert dataset_name == ll_term.pdf_base_config["name"], "Likelihood names do not match."
+                ll_term.set_data(d)
 
         self._data = data
 
@@ -79,13 +82,12 @@ def get_expectation_values(self, **kwargs) -> dict:
         Return total expectation values (summed over all likelihood terms with the same name)
         given a number of named parameters (kwargs)
         """
+        generate_values = self.parameters(**kwargs)  # kwarg or nominal value
         ret = dict()
-        # ancillary likelihood does not contribute
-        for ll in self._likelihood.likelihood_list[:-1]:
-
-            ll_pars = list(ll.rate_parameters.keys()) + list(ll.shape_parameters.keys())
-            ll_pars += ["livetime_days"]
-            call_args = {k: i for k, i in kwargs.items() if k in ll_pars}
+        # TODO: Make a self.likelihood_temrs dict with the likelihood names as keys and the corresponding likelihood terms as values.
+        for ll, parameter_names in zip(self._likelihood.likelihood_list[:-1],
+                                       self._likelihood.likelihood_parameters):  # ancillary likelihood does not contribute
+            call_args = {k: i for k, i in generate_values.items() if k in parameter_names}  # WARNING: This silently drops parameters it can't handle!
 
             mus = ll(full_output=True, **call_args)[1]
             for n, mu in zip(ll.source_name_list, mus):
@@ -117,7 +119,10 @@ def _build_ll_from_config(self, likelihood_config: dict) -> "LogLikelihoodSum":
             blueice_config["livetime_days"] = self.parameters[
                 blueice_config["livetime_parameter"]].nominal_value
             for p in self.parameters:
-                blueice_config[p.name] = blueice_config.get(p.name, p.nominal_value)
+                # adding the nominal rate values will screw things up in blueice!
+                # So here we're just adding the nominal values of all other parameters
+                if p.ptype != "rate":
+                    blueice_config[p.name] = blueice_config.get(p.name, p.nominal_value)
 
             # add all parameters to extra_dont_hash for each source unless it is used:
             for i, source in enumerate(config["sources"]):
@@ -150,25 +155,28 @@ def _build_ll_from_config(self, likelihood_config: dict) -> "LogLikelihoodSum":
                         "Only rate multipliers that end on _rate_multiplier"
                         " are currently supported.")
 
-                # Set shape parameters
-                shape_parameters = [
-                    p for p in source["parameters"] if self.parameters[p].ptype == "shape"]
-                if shape_parameters:
-                    # TODO: Implement setting shape parameters
-                    raise NotImplementedError("Shape parameters are not yet supported.")
-
                 # Set efficiency parameters
                 if source.get("apply_efficiency", False):
                     self._set_efficiency(source, ll)
 
+                # Set shape parameters
+                shape_parameters = [
+                    p for p in source["parameters"] if self.parameters[p].ptype == "shape"]
+                for p in shape_parameters:
+                    anchors = self.parameters[p].blueice_anchors
+                    if anchors is None:
+                        raise ValueError(f"Shape parameter {p} does not have any anchors.")
+                    ll.add_shape_parameter(p, anchors=anchors, log_prior=None)
+
             ll.prepare()
             lls.append(ll)
 
         # Ancillary likelihood
         ll = CustomAncillaryLikelihood(self.parameters.with_uncertainty)
         lls.append(ll)
 
-        return LogLikelihoodSum(lls, likelihood_weights=likelihood_config["likelihood_weights"])
+        likelihood_weights = likelihood_config.get("likelihood_weights", None)
+        return LogLikelihoodSum(lls, likelihood_weights=likelihood_weights)
 
     def _build_data_generators(self) -> list:
         # last one is AncillaryLikelihood
@@ -179,20 +187,21 @@ def _build_data_generators(self) -> list:
     def _ll(self, **generate_values) -> float:
         return self._likelihood(**generate_values)
 
-    def _generate_data(self, **generate_values) -> list:
+    def _generate_data(self, **generate_values) -> dict:
         # generate_values are already filtered and filled by the nominal values\
         # through the generate_data method in the parent class
-        science_data = self._generate_science_data(**generate_values)
+        data = self._generate_science_data(**generate_values)
         ancillary_keys = self.parameters.with_uncertainty.names
         generate_values_anc = {k: v for k, v in generate_values.items() if k in ancillary_keys}
-        ancillary_measurements = self._generate_ancillary_measurements(
+        data["ancillary_likelihood"] = self._generate_ancillary_measurements(
             **generate_values_anc)
-        return science_data + [ancillary_measurements] + [dict_to_structured_array(generate_values)]
+        data["generate_values"] = dict_to_structured_array(generate_values)
+        return data
 
-    def _generate_science_data(self, **generate_values) -> list:
-        science_data = [
-            gen.simulate(**generate_values) for gen in self.data_generators]
-        return science_data
+    def _generate_science_data(self, **generate_values) -> dict:
+        science_data = [gen.simulate(**generate_values)
+                        for gen in self.data_generators]
+        return dict(zip(self.likelihood_names, science_data))
 
     def _generate_ancillary_measurements(self, **generate_values) -> dict:
         ancillary_measurements = {}
@@ -262,10 +271,10 @@ def __init__(self, parameters: Parameters):
         parameter_list = self.parameters.names
 
         self.constraint_functions = self._get_constraint_functions()
-        super().__init__(
-            func=self.ancillary_likelihood_sum,
-            parameter_list=parameter_list,
-            config=self.parameters.nominal_values)
+        super().__init__(func=self.ancillary_likelihood_sum,
+                         parameter_list=parameter_list,
+                         config=self.parameters.nominal_values)
+        self.pdf_base_config["name"] = "ancillary_likelihood"
 
     @property
     def constraint_terms(self) -> dict:
@@ -283,10 +292,11 @@ def set_data(self, d: np.array):
             d (np.array): Data of ancillary measurements, stored as numpy array
         """
         # This results in shifted constraint terms.
-        if set(d.keys()) != set(self.parameters.names):
+        d_dict = structured_array_to_dict(d)
+        if set(d_dict.keys()) != set(self.parameters.names):
             raise ValueError(
                 "The data dict must contain all parameters as keys in CustomAncillaryLikelihood.")
-        self.constraint_functions = self._get_constraint_functions(**d)
+        self.constraint_functions = self._get_constraint_functions(**d_dict)
 
     def ancillary_likelihood_sum(self, evaluate_at: dict) -> float:
         """Return the sum of all constraint terms.

alea/simulators.py

@@ -57,7 +57,7 @@ def __init__(self, ll_term):
             source_histograms.append(mh.Histdd(bins=bins))
 
         self.ll = ll
-        logging.info("initing simulator, binned: " + str(binned))
+        logging.debug("initing simulator, binned: " + str(binned))
         self.binned = binned
         self.bincs = bincs
         self.direction_names = direction_names

alea/template_source.py

@@ -56,7 +56,7 @@ def build_histogram(self):
                 slice_axis_limits = sa.get(
                     "slice_axis_limits", [bes[0], bes[-1]])
                 if sum_axis:
-                    logging.info(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,
@@ -65,7 +65,7 @@ def build_histogram(self):
                     h.axis_names = axis_names
                 else:
                     logging.debug(f"Normalization before slicing: {h.n}.")
-                    logging.info(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,8 +126,8 @@ def build_histogram(self):
         # Shouldn't be in HistogramSource... anyway
         self._n_events_histogram = h.similar_blank_histogram()
 
-        h *= self.config.get('histogram_multiplier', 1)
-        logging.debug(f"Multiplying histogram with histogram_multiplier {self.config.get('histogram_multiplier', 1)}. Histogram is now normalised to {h.n}.")
+        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}.")
 
         # Convert h to density...
         if self.config.get('in_events_per_bin'):