TF error bands, bug fixes, minor improvements

src/rhalphalib/function.py

@@ -38,6 +38,60 @@ def matrix_poly(n: int):
     return np.identity(n + 1)
 
 
+def params_from_roofit(fitresult, param_names=None):
+    install_roofit_helpers()
+    names = [p.GetName() for p in fitresult.floatParsFinal()]
+    means = fitresult.valueArray()
+    cov = fitresult.covarianceArray()
+    if param_names is not None:
+        pidx = np.array([names.index(pname) for pname in param_names])
+        means = means[pidx]
+        cov = cov[np.ix_(pidx, pidx)]
+    return means, cov
+
+
+def sum_terms(params, coefs, shape):
+    return (params * coefs).sum(axis=1).reshape(shape)
+
+
+def compute_band(params, coefs, shape, cov):
+    """
+    This computation follows the linearized error propagation method used in RooAbsReal::plotOnWithErrorBand() and RooCurve::calcBandInterval()
+    err(x) = F(x,a) C_ab F(x,b)
+    where F(x,a) = (f(x,a+da) - f(x,a-da))/2 (numerical partial derivative) and C_ab is the correlation matrix
+    """
+    # compute correlation matrix
+    std = np.sqrt(np.diag(cov))
+    corr = cov / std[:, None] / std[None, :]
+    # compute gradients
+    plus_vars = []
+    minus_vars = []
+    for i, par in enumerate(params):
+        val = par
+        err = np.sqrt(cov[i, i])
+        # temporarily vary param
+        params[i] = val + err
+        plus_vars.append(sum_terms(params, coefs, shape))
+        params[i] = val - err
+        minus_vars.append(sum_terms(params, coefs, shape))
+        # reset to central value
+        params[i] = val
+    # flatten
+    plus_vars = np.stack(plus_vars).reshape(len(params), -1)
+    minus_vars = np.stack(minus_vars).reshape(len(params), -1)
+    F = (plus_vars - minus_vars) / 2.0
+
+    # final values
+    def proj(F):
+        return F.T @ corr @ F
+
+    central = sum_terms(params, coefs, shape)
+    err2 = np.apply_along_axis(proj, 1, F.T).reshape(shape)
+    lo = central + np.sqrt(err2)
+    hi = central - np.sqrt(err2)
+    return lo, hi
+
+
 class BasisPoly:
     """
     Construct a multidimensional polynomial
@@ -112,6 +166,8 @@ def __init__(
             else:
                 paramsq = DependentParameter("_".join([self.name] + ["%s_parsq%d" % (d, i) for d, i in zip(self._dim_names, ipar)]), "{0}*{0}", param)
                 self._params[ipar] = paramsq
+        # covariance is not populated until after fit is performed
+        self._cov = None
 
     @property
     def name(self):
@@ -145,12 +201,16 @@ def parameters(self, newparams):
             if pnew.intermediate:
                 pnew.intermediate = False
         self._params = newparams
+        # covariance is invalidated whenever parameters are changed
+        self._cov = None
 
     def update_from_roofit(self, fit_result, from_deco=False):
         par_names = sorted([p for p in fit_result.floatParsFinal().contentsString().split(",") if self.name in p])
-        par_results = {p: round(fit_result.floatParsFinal().find(p).getVal(), 3) for p in par_names}
+        means, cov = params_from_roofit(fit_result, par_names)
+        par_results = {p: round(means[i], 3) for i, p in enumerate(par_names)}
         for par in self._params.reshape(-1):
             par.value = par_results[par.name]
+        self._cov = cov
 
     def set_parvalues(self, parvalues):
         for par, new_val in zip(self._params.reshape(-1), parvalues):
@@ -167,12 +227,13 @@ def coefficients(self, *xvals):
             bpolyval = self._transform(bpolyval)
         return bpolyval
 
-    def __call__(self, *vals, nominal: bool = False):
+    def __call__(self, *vals, nominal: bool = False, errorband: bool = False):
         """Evaluate the polynomial at the given values
 
         Parameters:
             vals: a ndarray for each dimension's values to evaluate the polynomial at
             nominal: set true to evaluate nominal polynomial (rather than create DependentParameter objects)
+            errorband: set true to output error band along with nominal
         """
         if len(vals) != len(self._order):
             raise ValueError("Not all dimension values specified")
@@ -193,7 +254,14 @@ def __call__(self, *vals, nominal: bool = False):
         coefficients = self.coefficients(*xvals).reshape(-1, parameters.size)
         if nominal:
             parameters = np.vectorize(lambda p: p.value)(parameters)
-            return (parameters * coefficients).sum(axis=1).reshape(shape)
+            nominal_vals = sum_terms(parameters, coefficients, shape)
+            if errorband:
+                if self._cov is None:
+                    raise RuntimeError("Can only compute error band after covariance matrix loaded using update_from_roofit()")
+                band = compute_band(parameters, coefficients, shape, self._cov)
+                return nominal_vals, band
+            else:
+                return nominal_vals
 
         out = np.full(coefficients.shape[0], None)
         for i in range(coefficients.shape[0]):
@@ -241,12 +309,12 @@ def __init__(self, prefix: str, param_in: np.ndarray, param_cov: np.ndarray):
 
         _, s, v = np.linalg.svd(param_cov)
         self._transform = np.sqrt(s)[:, None] * v
-        self._parameters = np.array([NuisanceParameter(prefix + str(i), "param") for i in range(param_in.size)])
+        self._parameters = np.array([NuisanceParameter(prefix + str(i + 1), "param") for i in range(param_in.size)])
         self._correlated = np.full(self._parameters.shape, None)
         for i in range(self._parameters.size):
             coef = self._transform[:, i]
             order = np.argsort(np.abs(coef))
-            self._correlated[i] = np.sum(self._parameters[order] * coef[order]) + param_in[i]
+            self._correlated[i] = np.sum(coef[order] * self._parameters[order]) + param_in[i]
 
     @classmethod
     def fromRooFitResult(cls, prefix: str, fitresult, param_names: Optional[List[str]] = None):
@@ -258,14 +326,7 @@ def fromRooFitResult(cls, prefix: str, fitresult, param_names: Optional[List[str
             param_names: optional list of parameter names to include in the vector. If None,
                 all parameters in the fit result will be included.
         """
-        install_roofit_helpers()
-        names = [p.GetName() for p in fitresult.floatParsFinal()]
-        means = fitresult.valueArray()
-        cov = fitresult.covarianceArray()
-        if param_names is not None:
-            pidx = np.array([names.index(pname) for pname in param_names])
-            means = means[pidx]
-            cov = cov[np.ix_(pidx, pidx)]
+        means, cov = params_from_roofit(fitresult, param_names)
         out = cls(prefix, means, cov)
         if param_names is not None:
             for p, name in zip(out.correlated_params, param_names):
@@ -279,3 +340,7 @@ def parameters(self):
     @property
     def correlated_params(self):
         return self._correlated
+
+    @property
+    def correlated_str(self):
+        return "\n".join(p.formula(rendering=True) for p in self._correlated)

src/rhalphalib/model.py

@@ -210,7 +210,7 @@ def addSample(self, sample: Sample):
             raise ValueError("Naming convention requires beginning of sample %r name to be %s" % (sample, self.name))
         if self._observable is not None:
             if not sample.observable == self._observable:
-                raise ValueError("Sample %r has an incompatible observable with channel %r" % (sample, self))
+                raise ValueError("Sample %r has an incompatible observable with channel %r:\n  %r\n  %r" % (sample, self, sample.observable, self._observable))
             sample.observable = self._observable
         else:
             self._observable = sample.observable
@@ -236,7 +236,7 @@ def setObservation(self, obs, read_sumw2: bool = False):
         observable = Observable(obs_name, binning)
         if self._observable is not None:
             if not observable == self._observable:
-                raise ValueError("Observation has an incompatible observable with channel %r" % self)
+                raise ValueError("Observation has an incompatible observable with channel %r:\n  %r\n  %r" % (self, observable, self._observable))
         else:
             self._observable = observable
         if read_sumw2:

src/rhalphalib/parameter.py

@@ -17,9 +17,10 @@ def __init__(self, name: str, value):
         self._intermediate = False
 
     def __repr__(self):
-        return "<%s (%s) instance at 0x%x>" % (
+        return "<%s (%s, %s) instance at 0x%x>" % (
             self.__class__.__name__,
             self._name,
+            self._value,
             id(self),
         )
 
@@ -356,6 +357,14 @@ def __init__(self, name: str, binning: Iterable[float]):
         super(Observable, self).__init__(name, np.nan)
         self._binning = np.array(binning)
 
+    def __repr__(self):
+        return "<%s (%s, %s) instance at 0x%x>" % (
+            self.__class__.__name__,
+            self._name,
+            self._binning,
+            id(self),
+        )
+
     def __eq__(self, other):
         if isinstance(other, Observable) and self._name == other._name and np.allclose(self._binning, other._binning):
             return True

src/rhalphalib/template_morph.py

@@ -11,7 +11,7 @@ class AffineMorphTemplate(object):
     hist: a numpy-histogram-like tuple of (sumw, edges, varname)
     """
 
-    def __init__(self, hist: tuple[np.ndarray, np.ndarray, str]):
+    def __init__(self, hist: Tuple[np.ndarray, np.ndarray, str]):
         self.sumw = hist[0]
         self.edges = hist[1]
         self.varname = hist[2]

src/rhalphalib/util.py

@@ -20,7 +20,7 @@ def _to_numpy(hinput, read_sumw2=False):
         if read_sumw2 and hinput[3].size != hinput[1].size - 1:
             raise ValueError("Sumw2 array and binning array are incompatible in tuple {}".format(hinput))
         return hinput
-    elif "<class 'ROOT.TH1" in str(type(hinput)):
+    elif "<class 'ROOT.TH1" in str(type(hinput)) or "<class cppyy.gbl.TH1" in str(type(hinput)):
         sumw = np.zeros(hinput.GetNbinsX())
         sumw2 = np.zeros(hinput.GetNbinsX())
         binning = np.zeros(sumw.size + 1)
@@ -61,12 +61,15 @@ def _to_numpy(hinput, read_sumw2=False):
             return (sumw, binning, name, sumw2)
         return (sumw, binning, name)
     else:
-        raise ValueError("Cannot understand template type of %r" % hinput)
+        raise ValueError("Cannot understand template type %r of %r" % (type(hinput), hinput))
 
 
 def _to_TH1(sumw, binning, name):
     import ROOT
 
+    # avoid creating fake error values when setting bin content manually
+    oldDefaultSumw2 = ROOT.TH1.GetDefaultSumw2()
+    ROOT.TH1.SetDefaultSumw2(False)
     h = ROOT.TH1D(name, "template;%s;Counts" % name, binning.size - 1, binning)
     if isinstance(sumw, tuple):
         for i, (w, w2) in enumerate(zip(sumw[0], sumw[1])):
@@ -75,6 +78,7 @@ def _to_TH1(sumw, binning, name):
     else:
         for i, w in enumerate(sumw):
             h.SetBinContent(i + 1, w)
+    ROOT.TH1.SetDefaultSumw2(oldDefaultSumw2)
     return h
 
 

tests/test_rhalphalib.py

@@ -69,7 +69,7 @@ def test_rhalphabet(tmpdir):
         failObs = failCh.getObservation()
         qcdparams = np.array([rl.IndependentParameter("qcdparam_ptbin%d_msdbin%d" % (ptbin, i), 0) for i in range(msd.nbins)])
         sigmascale = 10.0
-        scaledparams = failObs * (1 + sigmascale / np.maximum(1.0, np.sqrt(failObs))) ** qcdparams
+        scaledparams = failObs * (1 + 1.0 / np.maximum(1.0, np.sqrt(failObs))) ** (qcdparams * sigmascale)
         fail_qcd = rl.ParametericSample("ptbin%dfail_qcd" % ptbin, rl.Sample.BACKGROUND, msd, scaledparams)
         failCh.addSample(fail_qcd)
         pass_qcd = rl.TransferFactorSample("ptbin%dpass_qcd" % ptbin, rl.Sample.BACKGROUND, tf_MCtempl_params[ptbin, :], fail_qcd)
@@ -167,8 +167,8 @@ def test_rhalphabet(tmpdir):
             initial_qcd -= sample.getExpectation(nominal=True)
         if np.any(initial_qcd < 0.0):
             raise ValueError("initial_qcd negative for some bins..", initial_qcd)
-        sigmascale = 10  # to scale the deviation from initial
-        scaledparams = initial_qcd * (1 + sigmascale / np.maximum(1.0, np.sqrt(initial_qcd))) ** qcdparams
+        sigmascale = 10.0
+        scaledparams = initial_qcd * (1 + 1.0 / np.maximum(1.0, np.sqrt(initial_qcd))) ** (qcdparams * sigmascale)
         fail_qcd = rl.ParametericSample("ptbin%dfail_qcd" % ptbin, rl.Sample.BACKGROUND, msd, scaledparams)
         failCh.addSample(fail_qcd)
         pass_qcd = rl.TransferFactorSample("ptbin%dpass_qcd" % ptbin, rl.Sample.BACKGROUND, tf_params[ptbin, :], fail_qcd)