feat: allow getExpectation(eval=True) for samples with independent params; allow shapeU

src/rhalphalib/model.py

@@ -59,6 +59,16 @@ def channels(self):
     def parameters(self):
         return reduce(set.union, (c.parameters for c in self), set())
 
+    def getParameter(self, name: str) -> Optional[IndependentParameter]:
+        for p in self.parameters:
+            if p.name == name:
+                return p
+        raise KeyError("Parameter %r not found in channel %r" % (name, self))
+
+    def reset(self):
+        for p in self.parameters:
+            p.reset()
+
     def addChannel(self, channel: "Channel"):
         """Add a channel to the model.
 
@@ -283,6 +293,16 @@ def samples(self):
     def parameters(self):
         return reduce(set.union, (s.parameters for s in self), set())
 
+    def getParameter(self, name: str) -> Optional[IndependentParameter]:
+        for p in self.parameters:
+            if p.name == name:
+                return p
+        raise KeyError("Parameter %r not found in channel %r" % (name, self))
+
+    def reset(self):
+        for p in self.parameters:
+            p.reset()
+
     @property
     def observable(self) -> Observable:
         """

src/rhalphalib/parameter.py

@@ -13,6 +13,7 @@ class Parameter:
     def __init__(self, name: str, value):
         self._name = name
         self._value = value
+        self._initial_value = value
         self._hasPrior = False
         self._intermediate = False
 
@@ -36,6 +37,9 @@ def name(self, name):
     def value(self):
         return self._value
 
+    def reset(self):
+        self._value = self._initial_value
+
     @property
     def intermediate(self):
         """
@@ -331,7 +335,7 @@ def value(self) -> float:
         return eval(self.formula().format(**{p.name: p.value for p in self.getDependents(deep=True)}))
 
     def formula(self, rendering=False):
-        return "((((0.1875**x - 0.625)*x*x + 0.9375)*x + 0.5)*(x > -1)*(x < 1) + 1*(x >= 1))".replace("x", "{%s}" % self.original_name)
+        return "(((0.1875*x*x - 0.625)*x*x + 0.9375)*x + 0.5)*(x > -1)*(x < 1) + 1*(x >= 1)".replace("x", "{%s}" % self.original_name)
 
     def renderRoofit(self, workspace):
         import ROOT

src/rhalphalib/sample.py

@@ -125,6 +125,9 @@ class TemplateSample(Sample):
         force_positive: if True, negative values in the template will be set to 0
     """
 
+    EffectNonMatchingThreshold = 0.85  # Threshold for the number of bins that are filled both in nominal and effect_up/down shape
+    EffectMagnitudeThreshold = 0.5  # Threshold for the effect magnitude to trigger a warning
+
     def __init__(self, name: str, sampletype: int, template, force_positive: bool = False):
         super(TemplateSample, self).__init__(name, sampletype)
         sumw2 = None
@@ -217,61 +220,69 @@ def setParamEffect(self, param: Parameter, effect_up, effect_down=None, scale=No
             raise ValueError(f"Parameter '{param.name}' already exists in sample '{self.name}': {sorted([p.name for p in self.parameters])}")
         if isinstance(effect_up, np.ndarray):
             if len(effect_up) != self.observable.nbins:
-                raise ValueError("effect_up has the wrong number of bins (%d, expected %d)" % (len(effect_up), self.observable.nbins))
+                raise ValueError(f"effect_up has the wrong number of bins ({len(effect_up)}, expected {self.observable.nbins})")
         elif isinstance(effect_up, numbers.Number):
             if "shape" in param.combinePrior:
                 effect_up = np.full(self.observable.nbins, effect_up)
         else:
             effect_up, binning, _ = _to_numpy(effect_up)
             if not np.array_equal(binning, self.observable.binning):
                 raise ValueError("effect_up has incompatible binning with sample %r" % self)
-            zerobins = (self._nominal <= 0.0) | (effect_up <= 0.0)
+            zerobins_nominal = (self._nominal <= 0.0) | np.isclose(self._nominal, 0.0)
+            zerobins_effect = (effect_up <= 0.0) | np.isclose(effect_up, 0.0)
+            _fraction_non_matching = sum(~zerobins_nominal == ~zerobins_effect) / len(zerobins_nominal)
+            if _fraction_non_matching <= self.EffectNonMatchingThreshold and not nowarn:
+                logging.warning(f"effect_up ({param.name}, {self._name}) has filled bins not matching with the nominal ({_fraction_non_matching * 100}%) - vertical morphing will be wonky.")
+            zerobins = zerobins_nominal | zerobins_effect
             effect_up[zerobins] = 1.0
             effect_up[~zerobins] /= self._nominal[~zerobins]
         if np.sum(effect_up * self._nominal) <= 0:
-            logging.warning("effect_up ({}, {}) has magnitude less than 0, skipping".format(param.name, self._name))
+            logging.warning(f"effect_up ({param.name}, {self._name}) has magnitude less than 0, skipping")
             # raise error instead?
             return
         elif effect_down is None and np.all(effect_up == 1.0):
-            logging.warning("effect_up ({}, {}) = 1 and has no effect, skipping".format(param.name, self._name))
+            logging.warning(f"effect_up ({param.name}, {self._name}) = 1 and has no effect, skipping")
             # raise error instead?
             return
         _weighted_effect_magnitude = np.sum(abs(effect_up - 1) * self._nominal) / np.sum(self._nominal)
-        if "shape" in param.combinePrior and _weighted_effect_magnitude > 0.5 and not nowarn:
+        if "shape" in param.combinePrior and _weighted_effect_magnitude > self.EffectMagnitudeThreshold and not nowarn:
             logging.warning(
-                "effect_up ({}, {}) has magnitude greater than 50% ({:.2f}%), you might be passing absolute values instead of relative".format(
-                    param.name, self._name, _weighted_effect_magnitude * 100
-                )
+                f"effect_up ({param.name}, {self._name}) has magnitude greater than {round(100 * self.EffectMagnitudeThreshold)}% ({_weighted_effect_magnitude * 100:.2f}%), you might be passing absolute values instead of relative"
             )
         self._paramEffectsUp[param] = effect_up
 
         if effect_down is not None:
             if isinstance(effect_down, np.ndarray):
                 if len(effect_down) != self.observable.nbins:
-                    raise ValueError("effect_down has the wrong number of bins (%d, expected %d)" % (len(effect_down), self.observable.nbins))
+                    raise ValueError(f"effect_down has the wrong number of bins ({len(effect_down)}, expected {self.observable.nbins})")
             elif isinstance(effect_down, numbers.Number):
                 if "shape" in param.combinePrior:
                     effect_down = np.full(self.observable.nbins, effect_down)
             else:
                 effect_down, binning, _ = _to_numpy(effect_down)
                 if not np.array_equal(binning, self.observable.binning):
                     raise ValueError("effect_down has incompatible binning with sample %r" % self)
-                zerobins = (self._nominal <= 0.0) | (effect_down <= 0.0)
+                zerobins_nominal = (self._nominal <= 0.0) | np.isclose(self._nominal, 0.0)
+                zerobins_effect = (effect_down <= 0.0) | np.isclose(effect_down, 0.0)
+                _fraction_non_matching = sum(~zerobins_nominal == ~zerobins_effect) / len(zerobins_nominal)
+                if _fraction_non_matching <= self.EffectNonMatchingThreshold and not nowarn:
+                    logging.warning(f"effect_down ({param.name}, {self._name}) has filled bins not matching with the nominal ({_fraction_non_matching * 100}%) - vertical morphing will be wonky.")
+                zerobins = zerobins_nominal | zerobins_effect
                 effect_down[zerobins] = 1.0
                 effect_down[~zerobins] /= self._nominal[~zerobins]
                 if np.sum(effect_down * self._nominal) <= 0:
                     # TODO: warning? this can happen regularly
                     # we might even want some sort of threshold
                     return
                 elif np.all(effect_up == 1.0) and np.all(effect_down == 1.0):
+                    logging.warning(f"effect_down ({param.name}, {self._name}) = 1 and has no effect, skipping")
+                    del self._paramEffectsUp[param]
                     # some sort of threshold might be useful here as well
                     return
             _weighted_effect_magnitude = np.sum(abs(effect_down - 1) * self._nominal) / np.sum(self._nominal)
-            if "shape" in param.combinePrior and _weighted_effect_magnitude > 0.5 and not nowarn:
+            if "shape" in param.combinePrior and _weighted_effect_magnitude > self.EffectMagnitudeThreshold and not nowarn:
                 logging.warning(
-                    "effect_down ({}, {}) has magnitude greater than 50% ({:.2f}%), you might be passing absolute values instead of relative".format(
-                        param.name, self._name, _weighted_effect_magnitude * 100
-                    )
+                    f"effect_down ({param.name}, {self._name}) has magnitude greater than {round(100 * self.EffectMagnitudeThreshold)}% ({_weighted_effect_magnitude * 100:.2f}%), you might be passing absolute values instead of relative"
                 )
             self._paramEffectsDown[param] = effect_down
         else:
@@ -380,8 +391,10 @@ def getExpectation(self, nominal: bool = False, eval: bool = False):
                     continue
                 if param in self._paramEffectScales:
                     param_scaled = param * self._paramEffectScales[param]
+                    param_scale = self._paramEffectScales[param]
                 else:
                     param_scaled = param
+                    param_scale = 1.0
                 if isinstance(effect_up, DependentParameter):
                     out = out * effect_up
                 elif self._paramEffectsDown[param] is None:
@@ -396,13 +409,11 @@ def getExpectation(self, nominal: bool = False, eval: bool = False):
                         raise NotImplementedError("per-bin effects for other nuisance parameter types")
                 else:
                     effect_down = self.getParamEffect(param, up=False)
-                    smoothStep = SmoothStep(param_scaled)
+                    smoothStep = SmoothStep(param) * param_scale
                     if param.combinePrior == "shape":
                         combined_effect = smoothStep * (1 + (effect_up - 1) * param_scaled) + (1 - smoothStep) * (1 - (effect_down - 1) * param_scaled)
-                    elif param.combinePrior == "shapeN":
-                        combined_effect = smoothStep * (effect_up**param_scaled) + (1 - smoothStep) / (effect_down**param_scaled)
-                    elif param.combinePrior == "lnN":
-                        # TODO: ensure scalar effect
+                    elif param.combinePrior in ["shapeN", "shapeU", "lnN"]:
+                        # TODO: ensure scalar effect for lnN
                         combined_effect = smoothStep * (effect_up**param_scaled) + (1 - smoothStep) / (effect_down**param_scaled)
                     else:
                         raise NotImplementedError("per-bin effects for other nuisance parameter types")
@@ -446,10 +457,14 @@ def renderRoofit(self, workspace):
                     continue
                 name = self.name + "_" + param.name + "Up"
                 shape = nominal * effect_up
+                if np.sum(shape) < 0:
+                    raise RuntimeError(f"Sample '{self.name}' has negative shape for parameter '{param.name}Up' with value {shape}. Can't build workspace.")
                 rooTemplate = ROOT.RooDataHist(name, name, ROOT.RooArgList(rooObservable), _to_TH1(shape, self.observable.binning, self.observable.name))
                 workspace.add(rooTemplate)
                 name = self.name + "_" + param.name + "Down"
                 shape = nominal * self.getParamEffect(param, up=False)
+                if np.sum(shape) < 0:
+                    raise RuntimeError(f"Sample '{self.name}' has negative shape for parameter '{param.name}Down with value {shape}. Can't build workspace.")
                 rooTemplate = ROOT.RooDataHist(name, name, ROOT.RooArgList(rooObservable), _to_TH1(shape, self.observable.binning, self.observable.name))
                 workspace.add(rooTemplate)