Skip to content

Commit

Permalink
Merge branch 'develop'
Browse files Browse the repository at this point in the history
  • Loading branch information
@WolfgangWaltenberger
WolfgangWaltenberger committed Sep 6, 2023
2 parents c951196 + 10ad9e0 commit 8299df4
Show file tree
Hide file tree
Showing 6 changed files with 54 additions and 87 deletions.
6 changes: 6 additions & 0 deletions ReleaseNotes
View file
Original file line number Diff line number Diff line change
@@ -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
=======================================================

Expand Down
6 changes: 6 additions & 0 deletions docs/manual/source/ReleaseUpdate.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 <analysesCombination>`
* added smodels version to output of "txt" printer

New in Version 2.3.1:
^^^^^^^^^^^^^^^^^^^^^

Expand Down
4 changes: 3 additions & 1 deletion smodels/tools/printer.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
118 changes: 34 additions & 84 deletions smodels/tools/pyhfInterface.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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):
"""
Expand Down Expand Up @@ -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:
Expand All @@ -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
Expand Down
5 changes: 4 additions & 1 deletion smodels/tools/simplifiedLikelihoods.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
2 changes: 1 addition & 1 deletion smodels/version
View file
Original file line number Diff line number Diff line change
@@ -1 +1 @@
2.3.1
2.3.2

0 comments on commit 8299df4

Please sign in to comment.