Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Removes all hash for parameters not used for each source, and for all… #37

Merged
merged 10 commits into from
Jul 17, 2023
47 changes: 43 additions & 4 deletions alea/blueice_extended_model.py
Original file line number Diff line number Diff line change
@@ -1,4 +1,6 @@
from pydoc import locate # to lookup likelihood class
from typing import List

from alea.statistical_model import StatisticalModel
from alea.simulators import BlueiceDataGenerator
from alea.utils import adapt_likelihood_config_for_blueice
Expand Down Expand Up @@ -76,7 +78,13 @@ def get_expectation_values(self, **kwargs) -> dict:
"""
ret = dict()
for ll in self._likelihood.likelihood_list[:-1]: # ancillary likelihood does not contribute
mus = ll(full_output=True, **kwargs)[1]

kdund marked this conversation as resolved.
Show resolved Hide resolved
ll_pars = list(ll.rate_parameters.keys()) + list(ll.shape_parameters.keys())
kdund marked this conversation as resolved.
Show resolved Hide resolved
call_args = {k:i for k, i in kwargs.items() if k in ll_pars}
dachengx marked this conversation as resolved.
Show resolved Hide resolved
if "livetime_days" in kwargs:
kdund marked this conversation as resolved.
Show resolved Hide resolved
call_args["livetime_days"] = kwargs["livetime_days"]
dachengx marked this conversation as resolved.
Show resolved Hide resolved

mus = ll(full_output=True, **call_args)[1]
for n, mu in zip(ll.source_name_list, mus):
ret[n] = ret.get(n, 0) + mu
return ret
Expand Down Expand Up @@ -105,13 +113,27 @@ def _build_ll_from_config(self, likelihood_config: dict) -> "LogLikelihoodSum":
config, template_folder_list)
blueice_config["livetime_days"] = self.parameters[
blueice_config["livetime_parameter"]].nominal_value
ll = likelihood_object(blueice_config)
for p in self.parameters:
blueice_config[p.name] = blueice_config.get(p.name, p.nominal_value)
kdund marked this conversation as resolved.
Show resolved Hide resolved

# add all parameters to extra_dont_hash for each source unless it is used:
for i, source in enumerate(config["sources"]):
parameters_to_ignore: List[str] = [p.name for p in self.parameters if (p.ptype == "shape")
dachengx marked this conversation as resolved.
Show resolved Hide resolved
dachengx marked this conversation as resolved.
Show resolved Hide resolved
dachengx marked this conversation as resolved.
Show resolved Hide resolved
dachengx marked this conversation as resolved.
Show resolved Hide resolved
& (p.name not in source["parameters"])]
dachengx marked this conversation as resolved.
Show resolved Hide resolved
dachengx marked this conversation as resolved.
Show resolved Hide resolved
# no efficiency affects PDF:
parameters_to_ignore += [ p.name for p in self.parameters if (p.ptype == "efficiency")]
dachengx marked this conversation as resolved.
Show resolved Hide resolved
dachengx marked this conversation as resolved.
Show resolved Hide resolved
dachengx marked this conversation as resolved.
Show resolved Hide resolved
dachengx marked this conversation as resolved.
Show resolved Hide resolved
parameters_to_ignore += source.get("extra_dont_hash_settings", [])

# ignore all shape parameters known to this model not named specifically in the source:
dachengx marked this conversation as resolved.
Show resolved Hide resolved
blueice_config["sources"][i]["extra_dont_hash_settings"] = parameters_to_ignore

ll = likelihood_object(blueice_config) # noqa: E127
dachengx marked this conversation as resolved.
Show resolved Hide resolved

for source in config["sources"]:

# Set rate parameters
rate_parameters = [
p for p in source["parameters"] if self.parameters[p].type == "rate"]
p for p in source["parameters"] if self.parameters[p].ptype == "rate"]
dachengx marked this conversation as resolved.
Show resolved Hide resolved
dachengx marked this conversation as resolved.
Show resolved Hide resolved
if len(rate_parameters) != 1:
raise ValueError(
f"Source {source['name']} must have exactly one rate parameter.")
Expand All @@ -125,11 +147,28 @@ def _build_ll_from_config(self, likelihood_config: dict) -> "LogLikelihoodSum":

# Set shape parameters
shape_parameters = [
p for p in source["parameters"] if self.parameters[p].type == "shape"]
p for p in source["parameters"] if self.parameters[p].ptype == "shape"]
dachengx marked this conversation as resolved.
Show resolved Hide resolved
dachengx marked this conversation as resolved.
Show resolved Hide resolved
if shape_parameters:
# TODO: Implement setting shape parameters
raise NotImplementedError("Shape parameters are not yet supported.")


# Set efficiency parameters
dachengx marked this conversation as resolved.
Show resolved Hide resolved
if source.get("apply_efficiency", False):
assert "efficiency_name" in source, "Unspecified efficiency_name for source {:s}".format(source["name"])
dachengx marked this conversation as resolved.
Show resolved Hide resolved
dachengx marked this conversation as resolved.
Show resolved Hide resolved
efficiency_name = source["efficiency_name"]
assert efficiency_name in source["parameters"], "The efficiency_name for source {:s} is not in its parameter list".format(source["name"])
dachengx marked this conversation as resolved.
Show resolved Hide resolved
dachengx marked this conversation as resolved.
Show resolved Hide resolved
efficiency_parameter = self.parameters[efficiency_name]
assert efficiency_parameter.ptype == "efficiency", "The parameter {:s} must" \
dachengx marked this conversation as resolved.
Show resolved Hide resolved
dachengx marked this conversation as resolved.
Show resolved Hide resolved
" be an efficiency".format(efficiency_name)
dachengx marked this conversation as resolved.
Show resolved Hide resolved
dachengx marked this conversation as resolved.
Show resolved Hide resolved
limits = efficiency_parameter.fit_limits
assert 0 <= limits[0], 'Efficiency parameters including {:s} must be' \
dachengx marked this conversation as resolved.
Show resolved Hide resolved
dachengx marked this conversation as resolved.
Show resolved Hide resolved
dachengx marked this conversation as resolved.
Show resolved Hide resolved
' constrained to be nonnegative'.format(efficiency_name)
dachengx marked this conversation as resolved.
Show resolved Hide resolved
assert np.isfinite(limits[1]), 'Efficiency parameters including {:s} must be' \
dachengx marked this conversation as resolved.
Show resolved Hide resolved
dachengx marked this conversation as resolved.
Show resolved Hide resolved
' constrained to be finite'.format(efficiency_name)
dachengx marked this conversation as resolved.
Show resolved Hide resolved
ll.add_shape_parameter(efficiency_name, anchors=(limits[0], limits[1]))


kdund marked this conversation as resolved.
Show resolved Hide resolved
ll.prepare()
lls.append(ll)

Expand Down
22 changes: 18 additions & 4 deletions alea/examples/unbinned_wimp_statistical_model.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,18 @@ parameter_definition:
- null
fit_guess: 1.0

signal_efficiency:
nominal_value: 1.0
ptype: efficiency
uncertainty: 0.1
relative_uncertainty: true
fittable: true
fit_limits:
- 0
- 10.
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
Expand Down Expand Up @@ -76,9 +88,10 @@ likelihood_config:
parameters:
- wimp_rate_multiplier
- wimp_mass
- signal_efficiency
template_filename: wimp50gev_template.h5
apply_efficiency: False
efficiency_name: 'signal_eff' # TODO: Check
apply_efficiency: True
efficiency_name: signal_efficiency # TODO: Check

# SR1
- name: sr1
Expand All @@ -104,6 +117,7 @@ likelihood_config:
parameters:
- wimp_rate_multiplier
- wimp_mass
- signal_efficiency
template_filename: wimp50gev_template.h5
apply_efficiency: False
efficiency_name: 'wimp_eff' # TODO: Check
apply_efficiency: True
efficiency_name: signal_efficiency # TODO: Check
4 changes: 2 additions & 2 deletions alea/parameters.py
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@ class Parameter:
name (str): The name of the parameter.
nominal_value (float, optional): The nominal value of the parameter.
fittable (bool, optional): Indicates if the parameter is fittable or always fixed.
ptype (str, optional): The type of the parameter.
ptype (str, optional): The ptype of the parameter.
uncertainty (float or str, optional):
The uncertainty of the parameter. If a string,
it can be evaluated as a numpy or scipy function to define non-gaussian constraints.
Expand Down Expand Up @@ -39,7 +39,7 @@ def __init__(
self.name = name
self.nominal_value = nominal_value
self.fittable = fittable
self.type = ptype
self.ptype = ptype
self._uncertainty = uncertainty
self.relative_uncertainty = relative_uncertainty
self.blueice_anchors = blueice_anchors
Expand Down
7 changes: 4 additions & 3 deletions alea/statistical_model.py
Original file line number Diff line number Diff line change
Expand Up @@ -163,7 +163,8 @@ def store_data(
kw = {'metadata': metadata} if metadata is not None else dict()
toydata_to_file(file_name, data_list, data_name_list, **kw)

def get_expectation_values(self):
def get_expectation_values(self, **parameter_values):
dachengx marked this conversation as resolved.
Show resolved Hide resolved
mus = self._ll(full_output=True, **parameter_values)
dachengx marked this conversation as resolved.
Show resolved Hide resolved
return NotImplementedError("get_expectation_values is optional to implement")
kdund marked this conversation as resolved.
Show resolved Hide resolved

@property
Expand Down Expand Up @@ -284,9 +285,9 @@ def _confidence_interval_checks(
"You must set parameter_interval_bounds in the parameter config"
" or when calling confidence_interval")

if parameter_of_interest.type == "rate":
if parameter_of_interest.ptype == "rate":
try:
if parameter_of_interest.type == "rate" and poi_name.endswith("_rate_multiplier"):
if parameter_of_interest.ptype == "rate" and poi_name.endswith("_rate_multiplier"):
source_name = poi_name.replace("_rate_multiplier", "")
else:
source_name = poi_name
Expand Down
Loading