From d91935e05461d56453e9fd40a254f1c0b8edc566 Mon Sep 17 00:00:00 2001 From: Wolfgang Waltenberger Date: Wed, 23 Aug 2023 10:27:54 +0200 Subject: [PATCH 1/4] use hessian also for 1d for computing sigma_mu --- smodels/tools/simplifiedLikelihoods.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/smodels/tools/simplifiedLikelihoods.py b/smodels/tools/simplifiedLikelihoods.py index 88187c7f6..fb086eb55 100755 --- a/smodels/tools/simplifiedLikelihoods.py +++ b/smodels/tools/simplifiedLikelihoods.py @@ -781,8 +781,11 @@ def lmax(self, return_nll=False, allowNegativeSignals=False): muhat = float(dn[0]) if abs(self.model.nsignal) > 1e-100: muhat = float(dn[0] / self.model.nsignal[0]) - sigma_mu = np.sqrt(self.model.observed[0] + self.model.covariance[0][0]) + #sigma_mu2 = np.sqrt(self.model.observed[0] / self.model.nsignal[0] + self.model.covariance[0][0] ) + theta_hat = self.findThetaHat( muhat ) + sigma_mu = self.getSigmaMu ( muhat, theta_hat[0] ) ret= self.likelihood( return_nll=return_nll, mu = muhat ) + # print ( "sigma_mu", sigma_mu, "old", sigma_mu2 ) return { "lmax": ret, "muhat": muhat, "sigma_mu": sigma_mu } fmh = self.findMuHat( allowNegativeSignals=allowNegativeSignals, extended_output=True, return_nll=return_nll From e3244e4002d8cb01adcf2712ea2eb88b2f03249f Mon Sep 17 00:00:00 2001 From: Wolfgang Waltenberger Date: Fri, 18 Aug 2023 17:27:12 +0200 Subject: [PATCH 2/4] three attempts at getting sigma_mu --- smodels/tools/pyhfInterface.py | 118 ++++++++++----------------------- 1 file changed, 34 insertions(+), 84 deletions(-) diff --git a/smodels/tools/pyhfInterface.py b/smodels/tools/pyhfInterface.py index a8cda5c65..267b25d31 100755 --- a/smodels/tools/pyhfInterface.py +++ b/smodels/tools/pyhfInterface.py @@ -592,76 +592,6 @@ def exponentiateNLL(self, twice_nll, doIt): return np.exp(-twice_nll / 2.0) return twice_nll / 2.0 - def getSigmaMu(self, workspace): - """given a workspace, compute a rough estimate of sigma_mu, - the uncertainty of mu_hat""" - obss, bgs, bgVars, nsig = {}, {}, {}, {} - channels = workspace.channels - for chdata in workspace["channels"]: - if not chdata["name"] in channels: - continue - bg = 0.0 - var = 0.0 - ns = 0.0 - for sample in chdata["samples"]: - if sample["name"] == "bsm": - ns = sample["data"][0] - else: - tbg = sample["data"][0] - bg += tbg - mult_factor_hi = 1. - mult_factor_lo = 1. - add_factor_hi = 0. - add_factor_lo = 0. - # To be tested - Don't take aux data into account - for modifier in sample["modifiers"]: - if modifier["type"] in ["normfactor","shapesys","staterror","shapefactor"]: - if type(modifier["data"]) == list: - mult_factor_hi *= modifier["data"][0] - mult_factor_lo *= modifier["data"][0] - elif modifier["type"] == "normsys": - if type(modifier["data"]) == dict: - mult_factor_hi *= modifier["data"]["hi"] - mult_factor_lo *= modifier["data"]["lo"] - elif modifier["type"] == "histosys": - # If many histosys modifiers, only the last one counts - if type(modifier["data"]) == dict: - add_factor_hi = modifier["data"]["hi_data"][0] - tbg - add_factor_lo = modifier["data"]["lo_data"][0] - tbg - elif modifier["type"] == "lumi": - for param in workspace["measurements"][0]["config"]["parameters"]: - if param["name"] == "lumi": - mult_factor_hi *= max(param["bounds"][0]) - mult_factor_lo *= min(param["bounds"][0]) - hi = mult_factor_hi*(tbg+add_factor_hi) - lo = mult_factor_lo*(tbg+add_factor_lo) - delta = max((hi - tbg, tbg - lo)) - var += delta**2 - nsig[chdata["name"]] = ns - bgs[chdata["name"]] = bg - bgVars[chdata["name"]] = var - for chdata in workspace["observations"]: - if not chdata["name"] in channels: - continue - obss[chdata["name"]] = chdata["data"][0] - vars = [] - for c in channels: - # poissonian error - if nsig[c]==0.: - nsig[c]=1e-5 - poiss = abs(obss[c]-bgs[c]) / nsig[c] - gauss = bgVars[c] / nsig[c]**2 - vars.append ( poiss + gauss ) - var_mu = np.sum ( vars ) - n = len ( obss ) - # print ( f" sigma_mu from pyhf uncorr {var_mu} {n} " ) - sigma_mu = float ( np.sqrt ( var_mu / (n**2) ) ) - # self.sigma_mu = sigma_mu - return sigma_mu - #import IPython - #IPython.embed() - #sys.exit() - def lmax( self, workspace_index=None, return_nll=False, expected=False, allowNegativeSignals=False): """ @@ -697,27 +627,46 @@ def lmax( self, workspace_index=None, return_nll=False, expected=False, # Same modifiers_settings as those used when running the 'pyhf cls' command line msettings = {"normsys": {"interpcode": "code4"}, "histosys": {"interpcode": "code4p"}} model = workspace.model(modifier_settings=msettings) + muhat, maxNllh = model.config.suggested_init(), float("nan") + sigma_mu = float("nan") # obs = workspace.data(model) try: bounds = model.config.suggested_bounds() if allowNegativeSignals: bounds[model.config.poi_index] = (-5., 10. ) - muhat, maxNllh = pyhf.infer.mle.fit(workspace.data(model), model, - return_fitted_val=True, par_bounds = bounds ) - if False: # get sigma_mu from hessian - pyhf.set_backend(pyhf.tensorlib, 'minuit') - muhat, maxNllh,o = pyhf.infer.mle.fit(workspace.data(model), model, - return_fitted_val=True, par_bounds = bounds, - return_result_obj = True ) - sigma_mu = float ( np.sqrt ( o.hess_inv[0][0] ) ) * self.scale - # print ( f"\n>>> sigma_mu from hessian {sigma_mu:.2f}" ) - pyhf.set_backend(pyhf.tensorlib, 'scipy') - - muhat = float ( muhat[model.config.poi_index]*self.scale ) + #import time + #t0 = time.time() + try: + muhat, maxNllh, o = pyhf.infer.mle.fit(workspace.data(model), model, + return_fitted_val=True, par_bounds = bounds, return_result_obj = True ) + sigma_mu = self.scale / float ( abs ( o.jac[model.config.poi_index] ) ) ## why did this not work? + except (pyhf.exceptions.FailedMinimization,ValueError) as e: + pass + #t1 = time.time() + #print ( f"first fit {t1-t0:.2f}s sigma_mu {sigma_mu:.3f}" ) + if hasattr ( o, "hess_inv" ): # maybe the backend gets changed + sigma_mu = float ( np.sqrt ( o.hess_inv[model.config.poi_index][model.config.poi_index] ) ) * self.scale + else: + sigma_mu_temp = 1. + try: + _1, _2, o = pyhf.infer.mle.fit(workspace.data(model), model, + return_fitted_val=True, return_result_obj = True, init_pars = list(muhat), method="BFGS" ) + sigma_mu_temp = float ( np.sqrt ( o.hess_inv[model.config.poi_index][model.config.poi_index] ) ) + except (pyhf.exceptions.FailedMinimization,ValueError) as e: + pass + if abs ( sigma_mu_temp - 1.0 ) > 1e-5: + sigma_mu = sigma_mu_temp * self.scale + else: + _, _, o = pyhf.infer.mle.fit(workspace.data(model), model, + return_fitted_val=True, return_result_obj = True, method="L-BFGS-B" ) + sigma_mu_temp = float ( np.sqrt ( o.hess_inv.todense()[model.config.poi_index][model.config.poi_index] ) ) + if abs ( sigma_mu_temp - 1.0 ) > 1e-8: + sigma_mu = sigma_mu_temp * self.scale +# sigma_mu = float ( o.unc[model.config.poi_index] ) * self.scale except (pyhf.exceptions.FailedMinimization, ValueError) as e: logger.error(f"pyhf mle.fit failed {e}") - muhat, maxNllh = float("nan"), float("nan") + muhat = float ( muhat[model.config.poi_index]*self.scale ) try: lmax = maxNllh.tolist() except: @@ -726,8 +675,9 @@ def lmax( self, workspace_index=None, return_nll=False, expected=False, lmax = float(lmax) except: lmax = float(lmax[0]) - sigma_mu = self.getSigmaMu ( workspace ) lmax = self.exponentiateNLL(lmax, not return_nll) + #t2 = time.time() + #print ( f"second fit {t2-t1:.2f}s sigma_mu {sigma_mu:.3f}" ) ret = { "lmax": lmax, "muhat": muhat, "sigma_mu": sigma_mu } self.data.cached_lmaxes[workspace_index] = ret return ret From 285a4eec64e1fef1e79db43f9d9833785dc0bb9d Mon Sep 17 00:00:00 2001 From: Wolfgang Waltenberger Date: Wed, 23 Aug 2023 11:40:34 +0200 Subject: [PATCH 3/4] preparing v2.3.2 --- ReleaseNotes | 6 ++++++ docs/manual/source/ReleaseUpdate.rst | 6 ++++++ smodels/tools/printer.py | 4 +++- 3 files changed, 15 insertions(+), 1 deletion(-) diff --git a/ReleaseNotes b/ReleaseNotes index 150feb573..d031b15ea 100644 --- a/ReleaseNotes +++ b/ReleaseNotes @@ -1,3 +1,9 @@ +Release v2.3.2, Wed 30 Aug 2023 +======================================================= + + * fixed bug in initialisation of analysis combination + * added smodels version to output of "txt" printer + Release v2.3.1, Mon 17 Jul 2023 ======================================================= diff --git a/docs/manual/source/ReleaseUpdate.rst b/docs/manual/source/ReleaseUpdate.rst index 77aa01203..79668bb0c 100644 --- a/docs/manual/source/ReleaseUpdate.rst +++ b/docs/manual/source/ReleaseUpdate.rst @@ -37,6 +37,12 @@ What's New ========== The major novelties of all releases since v1.0 are as follows: +New in Version 2.3.2: +^^^^^^^^^^^^^^^^^^^^^ + + * fixed bug in initialisation of :ref:`analyses combination ` + * added smodels version to output of "txt" printer + New in Version 2.3.1: ^^^^^^^^^^^^^^^^^^^^^ diff --git a/smodels/tools/printer.py b/smodels/tools/printer.py index bfa5b1111..8f7b9d8f9 100644 --- a/smodels/tools/printer.py +++ b/smodels/tools/printer.py @@ -333,8 +333,10 @@ def _formatOutputStatus(self, obj): for label in labels: par = obj.parameters[label] output += "# " + label + " = " + str(par) + '\n' + if obj.smodelsVersion: + output += f"# SModelS version: {obj.smodelsVersion}\n" if obj.databaseVersion: - output += "# Database version: %s\n" % obj.databaseVersion + output += f"# Database version: {obj.databaseVersion}\n" output += "=" * 80 + "\n" return output From 10ad9e06d60bb106adbe6ecffd8ff00e1cabaded Mon Sep 17 00:00:00 2001 From: Wolfgang Waltenberger Date: Wed, 23 Aug 2023 11:55:59 +0200 Subject: [PATCH 4/4] v232 --- smodels/version | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/smodels/version b/smodels/version index 2bf1c1ccf..f90b1afc0 100644 --- a/smodels/version +++ b/smodels/version @@ -1 +1 @@ -2.3.1 +2.3.2