First runner manipulating statistical model

alea/__init__.py

@@ -4,6 +4,8 @@
 
 from .models import *
 
+from .runner import *
+
 from .utils import *
 
 from .parameters import *

alea/examples/configs/unbinned_wimp_running.yaml

@@ -1,3 +1,4 @@
+statistical_model: alea.models.BlueiceExtendedModel
 statistical_model_config: unbinned_wimp_statistical_model.yaml
 
 poi: wimp_rate_multiplier
@@ -13,7 +14,7 @@ computation:
     parameters_in_common:
       {
         hypotheses: ["true", "null", "free"],
-        output_filename: "toymc_power_wimp_mass_{wimp_mass:d}_poi_expectation_{poi_expectation:.2f}.hdf5",
+        output_filename: "toymc_power_wimp_mass_{wimp_mass:d}_poi_expectation_{poi_expectation:.2f}.h5",
         n_mc: 5000,
         n_batch: 40,
       }
@@ -25,11 +26,11 @@ computation:
     parameters_in_common:
       {
         hypotheses: ["true", "null", "free"],
-        output_filename: "toymc_power_wimp_mass_{wimp_mass:d}_poi_expectation_{poi_expectation:.2f}.hdf5",
+        output_filename: "toymc_power_wimp_mass_{wimp_mass:d}_poi_expectation_{poi_expectation:.2f}.h5",
         n_mc: 5000,
         n_batch: 40,
       }
-    limit_threshold: "thresholds.hdf5"
+    limit_threshold: "thresholds.h5"
     toydata_mode: "generate"
     parameters_as_wildcards: ["poi_expectation", "n_mc", "n_batch"]
 
@@ -39,11 +40,11 @@ computation:
     parameters_in_common:
       {
         hypotheses: ["true", "null", "free"],
-        output_filename: "toymc_power_wimp_mass_{wimp_mass:d}_poi_expectation_{poi_expectation:.2f}.hdf5",
+        output_filename: "toymc_power_wimp_mass_{wimp_mass:d}_poi_expectation_{poi_expectation:.2f}.h5",
         n_mc: 5000,
         n_batch: 40,
         compute_confidence_interval: True,
-        limit_threshold: "thresholds.hdf5",
+        limit_threshold: "thresholds.h5",
       }
     toydata_mode: "generate"
 

alea/examples/gaussian_model.py

@@ -10,7 +10,7 @@ class GaussianModel(StatisticalModel):
     def __init__(
             self,
             parameter_definition: Optional[dict or list] = None,
-            **kwargs,
+            *args, **kwargs,
         ):
         """
         Initialise a model of a gaussian measurement (hatmu),
@@ -20,7 +20,7 @@ def __init__(
         """
         if parameter_definition is None:
             parameter_definition = ["mu", "sigma"]
-        super().__init__(parameter_definition=parameter_definition, **kwargs)
+        super().__init__(parameter_definition=parameter_definition, *args, **kwargs)
 
     def _ll(self, mu=None, sigma=None):
         hat_mu = self.data[0]['hat_mu'][0]

alea/model.py

@@ -58,6 +58,7 @@ def __init__(
             confidence_level: float = 0.9,
             confidence_interval_kind: str = "central",  # one of central, upper, lower
             confidence_interval_threshold: Callable[[float], float] = None,
+            **kwargs,
         ):
         """Initialize a statistical model"""
         if type(self) == StatisticalModel:
@@ -159,23 +160,36 @@ def data(self, data):
         self.is_data_set = True
 
     def store_data(
-            self, file_name, data_list, data_name_list=None, metadata = None):
+            self, file_name, data_list, data_name_list=None, metadata=None):
         """
-        Store a list of datasets (each on the form of a list of one or more structured arrays)
+        Store a list of datasets.
+        (each on the form of a list of one or more structured arrays or dicts)
         Using inference_interface, but included here to allow over-writing.
         structure would be: [[datasets1], [datasets2], ..., [datasetsn]]
         where each of datasets is a list of structured arrays
         if you specify, it is set, if not it will read from self.get_likelihood_term_names
         if not defined, it will be ["0", "1", ..., "n-1"]
         """
+        if all([isinstance(d, dict) for d in data_list]):
+            _data_list = [list(d.values()) for d in data_list]
+        elif all([isinstance(d, list) for d in data_list]):
+            _data_list = data_list
+        else:
+            raise ValueError(
+                'Unsupported mixed toydata format! '
+                'toydata should be a list of dict or a list of list',)
+
         if data_name_list is None:
             if hasattr(self, "likelihood_names"):
                 data_name_list = self.likelihood_names
             else:
-                data_name_list = ["{:d}".format(i) for i in range(len(data_list[0]))]
+                data_name_list = ["{:d}".format(i) for i in range(len(_data_list[0]))]
 
         kw = {'metadata': metadata} if metadata is not None else dict()
-        toydata_to_file(file_name, data_list, data_name_list, **kw)
+        if len(_data_list[0]) != len(data_name_list):
+            raise ValueError(
+                "The number of data sets and data names must be the same")
+        toydata_to_file(file_name, _data_list, data_name_list, **kw)
 
     def get_expectation_values(self, **parameter_values):
         return NotImplementedError("get_expectation_values is optional to implement")
@@ -195,7 +209,7 @@ def get_likelihood_term_from_name(self, likelihood_name):
         """
         if hasattr(self, "likelihood_names"):
             likelihood_names = self.likelihood_names
-            return {n:i for i,n in enumerate(likelihood_names)}[likelihood_name]
+            return {n: i for i, n in enumerate(likelihood_names)}[likelihood_name]
         else:
             raise NotImplementedError("The attribute likelihood_names is not defined.")
 

alea/models/blueice_extended_model.py

@@ -26,14 +26,17 @@ class BlueiceExtendedModel(StatisticalModel):
         likelihood_config (dict): A dictionary defining the likelihood.
     """
 
-    def __init__(self, parameter_definition: dict, likelihood_config: dict):
+    def __init__(
+            self,
+            parameter_definition: dict, likelihood_config: dict,
+            *args, **kwargs):
         """Initializes the statistical model.
 
         Args:
             parameter_definition (dict): A dictionary defining the model parameters.
             likelihood_config (dict): A dictionary defining the likelihood.
         """
-        super().__init__(parameter_definition=parameter_definition)
+        super().__init__(parameter_definition=parameter_definition, *args, **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")
@@ -62,14 +65,23 @@ def data(self) -> dict:
         return super().data
 
     @data.setter
-    def data(self, data: dict):
+    def data(self, data: dict or list):
         """
         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.
+
+        Args:
+            data (dict or list): Data of the statistical model.
+                If data is a list, it must be a list of length len(self.likelihood_names) + 1.
         """
         # iterate through all likelihood terms and set the science data in the blueice ll
         # last entry in data are the generate_values
+        if isinstance(data, list):
+            if len(data) != len(self.likelihood_names) + 1:
+                raise ValueError(
+                    f"Data must be a list of length {len(self.likelihood_names) + 1}")
+            data = dict(zip(self.likelihood_names + ["generate_values"], data))
         for i, (dataset_name, d) in enumerate(data.items()):
             if dataset_name != "generate_values":
                 ll_term = self._likelihood.likelihood_list[i]
@@ -138,7 +150,7 @@ def _build_ll_from_config(self, likelihood_config: dict) -> "LogLikelihoodSum":
             for i, source in enumerate(config["sources"]):
                 parameters_to_ignore: List[str] = [
                     p.name for p in self.parameters if (
-                        p.ptype == "shape") & (p.name not in source["parameters"])]
+                        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 += source.get("extra_dont_hash_settings", [])
@@ -158,6 +170,7 @@ def _build_ll_from_config(self, likelihood_config: dict) -> "LogLikelihoodSum":
                 rate_parameter = rate_parameters[0]
                 if rate_parameter.endswith("_rate_multiplier"):
                     rate_parameter = rate_parameter.replace("_rate_multiplier", "")
+                    # The ancillary term is handled in CustomAncillaryLikelihood
                     ll.add_rate_parameter(rate_parameter, log_prior=None)
                 else:
                     raise NotImplementedError(
@@ -175,6 +188,7 @@ def _build_ll_from_config(self, likelihood_config: dict) -> "LogLikelihoodSum":
                     anchors = self.parameters[p].blueice_anchors
                     if anchors is None:
                         raise ValueError(f"Shape parameter {p} does not have any anchors.")
+                    # The ancillary term is handled in CustomAncillaryLikelihood
                     ll.add_shape_parameter(p, anchors=anchors, log_prior=None)
 
             ll.prepare()
@@ -207,12 +221,23 @@ def _generate_data(self, **generate_values) -> dict:
         data["generate_values"] = dict_to_structured_array(generate_values)
         return data
 
+    def store_data(self, file_name, data_list, data_name_list=None, metadata=None):
+        """
+        Store data in a file.
+        Append the generate_values to the data_name_list.
+        """
+        if data_name_list is 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:
-        science_data = [gen.simulate(**generate_values)
-                        for gen in self.data_generators]
-        return dict(zip(self.likelihood_names, science_data))
+        """Generate the science data for all likelihood terms except the ancillary likelihood."""
+        science_data = [
+            gen.simulate(**generate_values) for gen in self.data_generators]
+        return dict(zip(self.likelihood_names[:-1], science_data))
 
     def _generate_ancillary_measurements(self, **generate_values) -> dict:
+        """Generate the ancillary measurements."""
         ancillary_measurements = {}
         anc_ll = self._likelihood.likelihood_list[-1]
         ancillary_generators = anc_ll._get_constraint_functions(**generate_values)

alea/parameters.py

@@ -1,5 +1,9 @@
 from typing import Any, Dict, List, Optional, Tuple
 
+# These imports are needed to evaluate the uncertainty string
+import numpy
+import scipy
+
 
 class Parameter:
     """
@@ -40,12 +44,12 @@ def __init__(
         self.nominal_value = nominal_value
         self.fittable = fittable
         self.ptype = ptype
-        self._uncertainty = uncertainty
+        self.uncertainty = uncertainty
         self.relative_uncertainty = relative_uncertainty
         self.blueice_anchors = blueice_anchors
         self.fit_limits = fit_limits
         self.parameter_interval_bounds = parameter_interval_bounds
-        self._fit_guess = fit_guess
+        self.fit_guess = fit_guess
         self.description = description
 
     def __repr__(self) -> str:
@@ -63,7 +67,7 @@ def uncertainty(self) -> float or Any:
         If the uncertainty is a string, it can be evaluated as a numpy or scipy function.
         """
         if isinstance(self._uncertainty, str):
-            # Evaluate the uncertainty if it's a string
+            # Evaluate the uncertainty if it's a string starting with "scipy." or "numpy."
             if self._uncertainty.startswith("scipy.") or self._uncertainty.startswith("numpy."):
                 return eval(self._uncertainty)
             else:
@@ -275,8 +279,8 @@ def __call__(
         if any(i is None for k, i in values.items()):
             emptypars = ", ".join([k for k, i in values.items() if i is None])
             raise AssertionError(
-                "All parameters must be set explicitly, or have a nominal value,"
-                " encountered for: " + emptypars)
+                "All parameters must be set explicitly, or have a nominal value, "
+                "not satisfied for: " + emptypars)
         return values
 
     def __getattr__(self, name: str) -> Parameter: