#!/usr/bin/env python3
"""
.. module:: analysesCombinations
:synopsis: a module with the methods required for computing the likelihood and
upper limits for a combination of theory predictions (analyses combination).
Mostly used by the TheoryPredictionsCombiner class.
.. moduleauthor:: Wolfgang Waltenberger <wolfgang.waltenberger@gmail.com>
.. moduleauthor:: Jamie Yellen <j.yellen.1@research.gla.ac.uk>
.. moduleauthor:: Andre Lessa <lessa.a.p@gmail.com>
"""
import numpy as np
from smodels.base.physicsUnits import fb, UnitXSec
from smodels.base.smodelsLogging import logger
from smodels.statistics.basicStats import CLsfromNLL, determineBrentBracket, \
findRoot, exponentiateNLL
import scipy.optimize as optimize
from smodels.statistics.exceptions import SModelSStatisticsError as SModelSError
from typing import Text, Tuple, Callable, Union, Dict
from smodels.tools.caching import roundCache, lru_cache
from smodels.statistics.basicStats import observed, apriori, aposteriori, \
NllEvalType
from smodels.matching.theoryPrediction import mu_digits
[docs]class AnaCombLikelihoodComputer(object):
def __init__(self, theoryPredictions: list, deltas_rel=None):
"""constructor.
:param theoryPredictions: the List of theory predictions
:param deltas_rel: relative uncertainty in signal (float). \
Default value is 20%.
"""
if len(theoryPredictions) == 0:
raise SModelSError("asking to combine zero theory predictions")
self.theoryPredictions = theoryPredictions
if deltas_rel is None:
from smodels.tools.runtime import _deltas_rel_default
deltas_rel = _deltas_rel_default
self.deltas_rel = deltas_rel
self.name = "anaComb"
[docs] @roundCache(argname='mu',argpos=1,digits=mu_digits)
def nll(
self,
mu: float = 1.0,
evaluationType : NllEvalType = observed,
asimov: Union[None,float] = None,
) -> float:
"""
Compute the negative log likelihood at a given mu. For now this is
just a proxy method to ease transition of codes that depend on SModelS.
:param mu: signal strength
:param evaluationType: one of: observed, apriori, aposteriori
:param asimov: if not None, compute llhd for asimov data with mu=asimov
"""
try:
mu = mu[0] # some of these methods use arrays with a single element
except:
pass
nll = 0.0
changed = False
for tp in self.theoryPredictions:
tmp = tp.nll(mu, evaluationType=evaluationType, asimov=asimov )
if tmp != None:
nll = nll + tmp #Add neg log llhds
changed = True
else:
return None
if changed == False:
return None
return nll
[docs] @lru_cache
def nll_min(
self,
allowNegativeSignals: bool = False,
evaluationType : NllEvalType = observed,
asimov: Union[None,float] = None,
) -> Union[Dict, None]:
"""find muhat and nll_min.
:param allowNegativeSignals: if true, then also allow for negative values
:param evaluationType: one of: observed, apriori, aposteriori
:returns: mu_hat, i.e. the maximum likelihood estimate of mu, if extended
output is requested, it returns a dictionary with mu_hat,
sigma_mu -- the standard deviation around mu_hat, and nll_min,
i.e. the likelihood at mu_hat
"""
muhats, weighted = [], []
totweight = 0.0
for tp in self.theoryPredictions:
muhat = asimov
if asimov is None:
muhat = tp.muhat(evaluationType=evaluationType)
sigma_mu = tp.sigma_mu(evaluationType=evaluationType)
if sigma_mu in [None, 0.0]:
sigma_mu = 1.0 # unity weights if no weights
if muhat != None:
muhats.append(muhat)
w = 1.0 / sigma_mu**2
weighted.append(w * muhat)
totweight += w
# for a single theory prediction, we return just that
if len(muhats)==1:
if muhat < 0. and not allowNegativeSignals:
muhat = 0.
retnllh = self.theoryPredictions[0].nll ( muhat, evaluationType = evaluationType, asimov=asimov )
ret = {"muhat": float(muhat), "sigma_mu": float(sigma_mu), "nll_min": retnllh }
return ret
if len(muhats) == 0:
logger.error(f"asked to compute muhat for combination, but no individual values")
ret = {"muhat": None, "sigma_mu": None, "nll_min": None }
return ret
toTry = [sum(weighted) / totweight]
#Add additional initialization points near the first one to make sure that minima is not lost due to small differences in the initialization values
fluc = max(1e-6, 1e-5*toTry[0]) #incase toTry[0] = 0.0, take 1e-06
toTry += [toTry[0] - fluc, toTry[0] + fluc]
def fun(mu):
# Make sure to always compute the correct llhd value (from theoryPrediction)
# and not used the cached value (which is constant for mu~=1 an mu~=0)
if isinstance(mu,np.ndarray):
x = float(mu[0])
else:
x = float(mu)
return self.nll(x, evaluationType=evaluationType, asimov=asimov)
if allowNegativeSignals:
toTry += [1.0, 0.0, 3.0, -1.0, 10.0, -3.0, 0.1, -0.1]
else:
toTry += [1.0, 0.0, 3.0, 10.0, 0.1]
for mu0 in toTry:
# Minimize with a delta_mu = 1e-3*mu0 when computing the first derivative
# (if delta_mu is too small, the derivative might give zero and terminate the
# minimization) o = scipy.optimize.minimize(fun, mu0 ) ## the unbounded method
mxm = float(max(muhats))
upper = 3.0 * mxm
if upper <= 1.5: #
upper = 1.5
bounds = [0, upper]
if allowNegativeSignals:
m = float(min(muhats))
if m < 0.0:
bounds[0] = 3.0 * m
bounds = [
tuple(bounds),
]
# if bounds[1] < bounds[0]:
# logger.error ( f"bounds are reversed, this should not happen" )
o = optimize.minimize(fun, mu0, bounds=bounds, tol=1e-9)
if not o.success:
logger.debug(
f"combiner.nll_min did not terminate successfully: {o.message} "
f"mu_hat={o.x} hess={o.hess_inv}"
)
# the inverted hessian is a good approximation for the variance at the
# minimum
invh = o.hess_inv
try:
invh = invh.todense()
except AttributeError as e:
pass
hessian = invh[0][0]
nll_ = o.fun
#removed if condition: "and nll_ < 998.0", because sometimes nll_ > 998.0. Not necessarily a bad thing?
#However lmax = np.exp(-nll_) for large nll -> 0
if hessian > 0.0 and o.success: # found a maximum. all good. #and nll_ < 998.0
break
# the hessian is negative meaning we found a maximum, not a minimum
if hessian <= 0.0:
logger.debug(
f"combiner.nll_min the hessian {hessian} is negative at mu_hat={o.x}. "+
"try again with different initialisation."
)
mu_hat = o.x[0]
# lmax = np.exp(-o.fun) # fun is *always* nll
if hessian < 0.0: #or nll_ > 998.0: again remove nll_ cut
logger.error(
"tried with several starting points to find maximum, always ended up in minimum. "
"bailing out."
)
return None
if not allowNegativeSignals and mu_hat < 0.0:
mu_hat = 0.0 # fixme for this case we should reevaluate the hessian!
sigma_mu = np.sqrt(hessian)
ret = {"muhat": float(mu_hat), "sigma_mu": float(sigma_mu) }
ret["nll_min"] = nll_
return ret
[docs] @lru_cache
def getUpperLimitOnMu(self, evaluationType : NllEvalType=observed,
allowNegativeSignals = False, nSigma : int = 0 ) -> float:
"""get upper limit on signal strength multiplier, i.e. value for mu for \
which CLs = 0.95
:param evaluationType: one of: observed, apriori, aposteriori
:param nSigma: the upper limit for central value (0),
+ 1 sigma, - 1 sigma, etc. For error bands.
:returns: upper limit on signal strength multiplier mu
"""
mu_hat, sigma_mu, clsRoot = self.getCLsRootFunc(
evaluationType=evaluationType,
allowNegativeSignals=allowNegativeSignals,
nSigma = nSigma )
a, b = determineBrentBracket(mu_hat, sigma_mu, clsRoot,
allowNegative = allowNegativeSignals )
mu_lim = findRoot(clsRoot, a, b, rtol=1e-03, xtol=1e-06 )
return mu_lim
[docs] def getTotalXSec ( self ):
xsec = 0.0*fb
for tp in self.theoryPredictions:
xsec += tp.xsection
return xsec
[docs] def getCLsRootFunc(self, evaluationType : NllEvalType=observed,
allowNegativeSignals : bool = False, nSigma : int = 0 ) \
-> Tuple[float, float, Callable]:
"""
Obtain the function "CLs-alpha[0.05]" whose root defines the upper limit,
plus mu_hat and sigma_mu
:param nSigma: the upper limit for central value (0),
+ 1 sigma, - 1 sigma, etc. For error bands.
:param evaluationType: one of: observed, apriori, aposteriori
"""
fmh = self.nll_min( evaluationType=evaluationType,
allowNegativeSignals=allowNegativeSignals )
mu_hat, sigma_mu = fmh["muhat"], fmh["sigma_mu"]
mu_hat = mu_hat if mu_hat is not None else 0.0
nll0 = self.nll(mu_hat, evaluationType=evaluationType,
asimov=None )
# a posteriori evaluationType is needed here
# mu_hat is mu_hat for signal_rel
fmh = self.nll_min(evaluationType=evaluationType,
allowNegativeSignals=allowNegativeSignals,
asimov = 0. )
nll0A = fmh["nll_min"]
# logger.error ( f"COMB nll0A {nll0A:.3f} mu_hatA {mu_hatA:.3f}" )
# return 1.
def clsRoot(mu: float, return_type: Text = "CLs-alpha") -> float:
# at - infinity this should be .95,
# at + infinity it should -.05
# Make sure to always compute the correct llhd value (from theoryPrediction)
# and not used the cached value (which is constant for mu~=1 an mu~=0)
nll = self.nll (mu,
evaluationType=evaluationType, asimov = None)
nllA = self.nll ( mu, evaluationType=evaluationType,
asimov = 0. )
if nll is None or nllA is None:
return None
ret = CLsfromNLL( nllA, nll0A, nll, nll0, (mu_hat>mu),
return_type=return_type, nSigma=nSigma )
return ret
return mu_hat, sigma_mu, clsRoot
[docs] @roundCache(argname='mu',argpos=1,digits=mu_digits)
def CLs( self, mu : float = 1., evaluationType : NllEvalType=observed,
return_type: Text = "CLs", nSigma : int = 0,
pmSigma : int = 0 ) -> float:
"""
Compute the exclusion confidence level of the model
:param mu: compute for the parameter of interest mu
:param evaluationType: one of: observed, apriori, aposteriori
:param return_type: (Text) can be "CLs-alpha", "1-CLs", "CLs" \
CLs-alpha: returns CLs - 0.05 \
alpha-CLs: returns 0.05 - CLs \
1-CLs: returns 1-CLs value \
CLs: returns CLs value
"""
assert pmSigma == 0, f"pmSigma {pmSigma} not supported by combiner"
assert return_type in ["CLs-alpha", "alpha-CLs", "1-CLs", "CLs"], \
f"Unknown return type: {return_type}."
_, _, clsRoot = self.getCLsRootFunc(evaluationType=evaluationType,
nSigma = nSigma )
return float(clsRoot(mu, return_type=return_type))
[docs] def getLlhds( self, muvals : list, evaluationType : NllEvalType=observed,
normalize : bool = True) -> dict:
"""
Compute the likelihoods for the individual analyses and the combined
likelihood.
Returns a dictionary with the analysis IDs as keys and the
likelihood values as values.
:param muvals: List with values for the signal strenth for which the
likelihoods must be evaluated.
:param evaluationType: one of: observed, apriori, aposteriori
:param normalize: If True normalizes the likelihood by its integral
over muvals.
"""
llhds = {}
llhds['combined'] = np.array(\
[exponentiateNLL(self.nll(mu,evaluationType=evaluationType)) for mu in muvals])
tpreds = self.theoryPredictions
for t in tpreds:
Id = t.analysisId()
t.computeStatistics( evaluationType = evaluationType )
l = np.array(\
[exponentiateNLL(t.nll(mu,evaluationType=evaluationType)) for mu in muvals])
llhds[Id]=l
if normalize:
# Compute delta mu
dmuvals = np.diff(muvals)
dmuvals = np.append(dmuvals,dmuvals[-1])
for Id,l in llhds.items():
# Compute norm (integral over mu)
norm = np.sum(l*dmuvals)
llhds[Id] = l/norm
return llhds