From 19fba6b343b6e65c464a3f1c4c6a4aadcb083242 Mon Sep 17 00:00:00 2001 From: Oliver Majersky Date: Sat, 9 Feb 2019 15:32:50 +0100 Subject: [PATCH 1/7] Add hacky way to specify object syst NP with uniform prior instead of Gaussian --- fbu/PyFBU.py | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/fbu/PyFBU.py b/fbu/PyFBU.py index 8cf0ac4..1a0aa23 100644 --- a/fbu/PyFBU.py +++ b/fbu/PyFBU.py @@ -28,6 +28,7 @@ def __init__(self,data=[],response=[],background={}, # [unfolding model parameters] self.prior = 'Uniform' self.priorparams = {} + self.obj_syst_flatprior = {'key': '__flat__', 'lower':-5, 'upper':5} self.regularization = regularization # [input] self.data = data # data list @@ -114,9 +115,11 @@ def run(self): bckgnuisances = mc.math.stack(bckgnuisances) if nobjsyst>0: - objnuisances = [ mc.Normal('gaus_%s'%name,mu=0.,tau=1.0#, - #observed=(True if self.systfixsigma!=0 else False) - ) + objnuisances = [ mc.Uniform('flat_%s'%name, + lower=self.obj_syst_flatprior['lower'], + upper=self.obj_syst_flatprior['upper']) + if self.obj_syst_flatprior['key'] in name else + mc.Normal('gaus_%s'%name,mu=0.,tau=1.0) for name in objsystkeys] objnuisances = mc.math.stack(objnuisances) @@ -158,7 +161,7 @@ def unfold(): print(self.nuts_kwargs) - + if self.mode: map_estimate = mc.find_MAP(model=model, method=self.MAP_method) print (map_estimate) @@ -190,7 +193,10 @@ def unfold(): #self.nuisancestrace[name] = copy.deepcopy(trace['gaus_%s'%name][:]) for name in objsystkeys: if self.systfixsigma==0.: - self.nuisancestrace[name] = trace['gaus_%s'%name][:] + if self.obj_syst_flatprior['key'] in name: + self.nuisancestrace[name] = trace['flat_%s'%name][:] + else: + self.nuisancestrace[name] = trace['gaus_%s'%name][:] #self.nuisancestrace[name] = copy.deepcopy(trace['gaus_%s'%name][:]) if self.monitoring: From 7a04f22e8523de3d475516e5af9772de77041467 Mon Sep 17 00:00:00 2001 From: Oliver Majersky Date: Sat, 15 Jun 2019 19:06:30 +0200 Subject: [PATCH 2/7] Add preliminary support to freeze NPs to a specific value --- fbu/PyFBU.py | 79 ++++++++++++++++++++++++++++++++++++++-------------- 1 file changed, 58 insertions(+), 21 deletions(-) diff --git a/fbu/PyFBU.py b/fbu/PyFBU.py index 1a0aa23..b2ae55a 100644 --- a/fbu/PyFBU.py +++ b/fbu/PyFBU.py @@ -29,6 +29,7 @@ def __init__(self,data=[],response=[],background={}, self.prior = 'Uniform' self.priorparams = {} self.obj_syst_flatprior = {'key': '__flat__', 'lower':-5, 'upper':5} + self.freeze_NPs = {} # nuisance parameters for which to fix value (not sampled) self.regularization = regularization # [input] self.data = data # data list @@ -94,6 +95,11 @@ def run(self): model = mc.Model() from .priors import wrapper + add_kwargs = dict() + if len(self.freeze_NPs) > 0: + print('Freezing values of following NPs:') + for key, val in self.freeze_NPs.items(): + print('{0}: {1}'.format(key, val)) with model: truth = wrapper(priorname=self.prior, low=self.lower,up=self.upper, @@ -102,25 +108,44 @@ def run(self): if nbckg>0: bckgnuisances = [] for name,err in zip(backgroundkeys,backgroundnormsysts): + try: + add_kwargs['observed'] = self.freeze_NPs[name] + except KeyError: + add_kwargs = dict() if err<0.: bckgnuisances.append( - mc.Uniform('norm_%s'%name,lower=0.,upper=3.) + mc.Uniform('norm_%s'%name,lower=0.,upper=3., **add_kwargs) ) else: - BoundedNormal = mc.Bound(mc.Normal, lower=(-1.0/err if err>0.0 else -inf)) - bckgnuisances.append( - BoundedNormal('gaus_%s'%name, - mu=0.,tau=1.0) + # for fixed NP, one cannot use observed in bounded + # distribution, so we have to use unbounded one + if 'observed' in add_kwargs: + bckgnuisances.append( + mc.Normal('gaus_%s'%name, mu=0.,tau=1.0, + **add_kwargs) + ) + else: + BoundedNormal = mc.Bound(mc.Normal, lower=(-1.0/err if err>0.0 else -inf)) + bckgnuisances.append( + BoundedNormal('gaus_%s'%name, mu=0.,tau=1.0) ) bckgnuisances = mc.math.stack(bckgnuisances) if nobjsyst>0: - objnuisances = [ mc.Uniform('flat_%s'%name, + objnuisances = list() + for name in objsystkeys: + try: + add_kwargs['observed'] = self.freeze_NPs[name] + except KeyError: + add_kwargs = dict() + if self.obj_syst_flatprior['key'] in name: + objnuisances.append(mc.Uniform('flat_%s'%name, lower=self.obj_syst_flatprior['lower'], - upper=self.obj_syst_flatprior['upper']) - if self.obj_syst_flatprior['key'] in name else - mc.Normal('gaus_%s'%name,mu=0.,tau=1.0) - for name in objsystkeys] + upper=self.obj_syst_flatprior['upper'], + **add_kwargs)) + else: + objnuisances.append(mc.Normal('gaus_%s'%name,mu=0., + tau=1.0, **add_kwargs)) objnuisances = mc.math.stack(objnuisances) # define potential to constrain truth spectrum @@ -185,19 +210,31 @@ def unfold(): self.nuisancestrace = {} if nbckg>0: for name,err in zip(backgroundkeys,backgroundnormsysts): - if err<0.: - self.nuisancestrace[name] = trace['norm_%s'%name][:] - #self.nuisancestrace[name] = copy.deepcopy(trace['norm_%s'%name][:]) - if err>0.: - self.nuisancestrace[name] = trace['gaus_%s'%name][:] - #self.nuisancestrace[name] = copy.deepcopy(trace['gaus_%s'%name][:]) + try: + if err<0.: + self.nuisancestrace[name] = trace['norm_%s'%name][:] + #self.nuisancestrace[name] = copy.deepcopy(trace['norm_%s'%name][:]) + if err>0.: + self.nuisancestrace[name] = trace['gaus_%s'%name][:] + #self.nuisancestrace[name] = copy.deepcopy(trace['gaus_%s'%name][:]) + except KeyError as e: + try: + tmp = self.freeze_NPs[name] + except KeyError as e: + print('Warning: Missing NP trace', e) for name in objsystkeys: if self.systfixsigma==0.: - if self.obj_syst_flatprior['key'] in name: - self.nuisancestrace[name] = trace['flat_%s'%name][:] - else: - self.nuisancestrace[name] = trace['gaus_%s'%name][:] - #self.nuisancestrace[name] = copy.deepcopy(trace['gaus_%s'%name][:]) + try: + if self.obj_syst_flatprior['key'] in name: + self.nuisancestrace[name] = trace['flat_%s'%name][:] + else: + self.nuisancestrace[name] = trace['gaus_%s'%name][:] + #self.nuisancestrace[name] = copy.deepcopy(trace['gaus_%s'%name][:]) + except KeyError as e: + try: + tmp = self.freeze_NPs[name] + except KeyError as e: + print('Warning: Missing NP trace', e) if self.monitoring: from fbu import monitoring From 9fa924d604295c3b4dea91e1fc61461cc7764301 Mon Sep 17 00:00:00 2001 From: Oliver Majersky Date: Mon, 8 Jul 2019 10:52:13 +0200 Subject: [PATCH 3/7] Add ability to smear backgrounds according to their MC stat --- fbu/PyFBU.py | 28 +++++++++++++++++++++++----- 1 file changed, 23 insertions(+), 5 deletions(-) diff --git a/fbu/PyFBU.py b/fbu/PyFBU.py index b2ae55a..222a9a2 100644 --- a/fbu/PyFBU.py +++ b/fbu/PyFBU.py @@ -36,8 +36,10 @@ def __init__(self,data=[],response=[],background={}, self.response = response # response matrix self.background = background # background dict self.backgroundsyst = backgroundsyst + self.backgrounderr = {} self.objsyst = objsyst self.systfixsigma = 0. + self.smear_bckgs = {} # backgrounds to be smeared in PE (according to MC stats) # [settings] self.rndseed = rndseed self.verbose = verbose @@ -59,14 +61,30 @@ def checklen(list1,list2): for bin in [self.lower,self.upper]: checklen(bin,responsetruthbins) #__________________________________________________________ - def fluctuate(self, data): + def fluctuate(self, data, err=None): random.seed(self.rndseed) - return random.poisson(data) + if err is None: + return random.poisson(data) + else: + return random.normal(data, err) #__________________________________________________________ def run(self): self.validateinput() - data = self.data - data = self.fluctuate(data) if self.rndseed>=0 else data + data = copy.deepcopy(self.data) + background = copy.deepcopy(self.background) + if len(self.smear_bckgs) > 0: + if self.rndseed >= 0: + for bckg in self.smear_bckgs: + try: + background[bckg] = self.fluctuate(background[bckg], + self.backgrounderr[bckg]) + except KeyError as e: + print('Error when trying to smear background {0}') + print('Check that the background exists in background' + ' and backgrounderr dictionaries.') + raise + else: + data = self.fluctuate(data) if self.rndseed>=0 else data # unpack background dictionaries backgroundkeys = self.backgroundsyst.keys() @@ -75,7 +93,7 @@ def run(self): backgrounds = [] backgroundnormsysts = array([]) if nbckg>0: - backgrounds = array([self.background[key] for key in backgroundkeys]) + backgrounds = array([background[key] for key in backgroundkeys]) backgroundnormsysts = array([self.backgroundsyst[key] for key in backgroundkeys]) # unpack object systematics dictionary From 6ab334ba352a641da2fa513a5fd1ebbf6952e89b Mon Sep 17 00:00:00 2001 From: Oliver Majersky Date: Sat, 20 Jul 2019 20:10:04 +0200 Subject: [PATCH 4/7] Preliminary support for gamma NPs for background MC stat uncertainty. --- fbu/PyFBU.py | 41 ++++++++++++++++++++++++++++++++++++++++- 1 file changed, 40 insertions(+), 1 deletion(-) diff --git a/fbu/PyFBU.py b/fbu/PyFBU.py index 222a9a2..754b849 100644 --- a/fbu/PyFBU.py +++ b/fbu/PyFBU.py @@ -1,5 +1,5 @@ import pymc3 as mc -from numpy import random, dot, array, inf +from numpy import random, dot, array, inf, sum, sqrt, reciprocal import theano import copy @@ -38,6 +38,9 @@ def __init__(self,data=[],response=[],background={}, self.backgroundsyst = backgroundsyst self.backgrounderr = {} self.objsyst = objsyst + self.include_gammas = False + self.gammas_treshold = 0.005 + self.nbins = 0 self.systfixsigma = 0. self.smear_bckgs = {} # backgrounds to be smeared in PE (according to MC stats) # [settings] @@ -60,6 +63,10 @@ def checklen(list1,list2): checklen(self.data,bin) for bin in [self.lower,self.upper]: checklen(bin,responsetruthbins) + + if self.include_gammas: + assert self.backgrounderr != {},\ + 'To include gammas, must provide background MC stat uncertainties' #__________________________________________________________ def fluctuate(self, data, err=None): random.seed(self.rndseed) @@ -89,13 +96,24 @@ def run(self): # unpack background dictionaries backgroundkeys = self.backgroundsyst.keys() nbckg = len(backgroundkeys) + self.nbins = len(background[next(iter(background))]) backgrounds = [] + backgrounds_err = [] backgroundnormsysts = array([]) if nbckg>0: backgrounds = array([background[key] for key in backgroundkeys]) + # backgrounds_err_sq = array([self.backgrounderr[key]**2 for key in backgroundkeys]) + backgrounds_err_sq = array([self.backgrounderr[key] for key in backgroundkeys]) + backgrounds_err_sq = backgrounds_err_sq**2 backgroundnormsysts = array([self.backgroundsyst[key] for key in backgroundkeys]) + # need summed total background and it's error for gamma NPs + # to take into account MC stat uncertainty of backgrounds + if self.include_gammas: + totalbckg = sum(backgrounds, axis=0) + totalbckg_err = sqrt(sum(backgrounds_err_sq, axis=0)) + # unpack object systematics dictionary objsystkeys = self.objsyst['signal'].keys() nobjsyst = len(objsystkeys) @@ -166,6 +184,20 @@ def run(self): tau=1.0, **add_kwargs)) objnuisances = mc.math.stack(objnuisances) + if self.include_gammas and nbckg > 0: + gammas = [mc.Uniform('flat_gamma_{0}'.format(i), lower=0., + upper=10.) for i in range(self.nbins)] + gammas = mc.math.stack(gammas) + + # construct the Poisson constraint on gammas + tau = (totalbckg/totalbckg_err)**2 + m = reciprocal(tau) + print(sqrt(m)) + print(reciprocal(sqrt(totalbckg))) + gamma_poissons = [ + mc.Poisson('poisson_gamma_{0}'.format(i), + mu=tau[i], observed=m[i]) for i in range(self.nbins)] + # define potential to constrain truth spectrum if self.regularization: truthpot = self.regularization.getpotential(truth) @@ -185,6 +217,9 @@ def unfold(): bckg = theano.dot(1. + bckgnuisances*bckgnormerr,smearedbackgrounds) + if self.include_gammas: + bckg = bckg * gammas + tresmat = array(resmat) reco = theano.dot(truth, tresmat) out = reco @@ -240,6 +275,10 @@ def unfold(): tmp = self.freeze_NPs[name] except KeyError as e: print('Warning: Missing NP trace', e) + if self.include_gammas: + for bin in range(self.nbins): + self.nuisancestrace['gamma_{0}'.format(bin)]\ + = trace['flat_gamma_{0}'.format(bin)][:] for name in objsystkeys: if self.systfixsigma==0.: try: From 4201300cb10a8dfff848c2043e9d66d676acb3d5 Mon Sep 17 00:00:00 2001 From: Oliver Majersky Date: Sat, 20 Jul 2019 23:12:57 +0200 Subject: [PATCH 5/7] Fix bug in Poisson constraint for gammas --- fbu/PyFBU.py | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/fbu/PyFBU.py b/fbu/PyFBU.py index 754b849..f2b69cd 100644 --- a/fbu/PyFBU.py +++ b/fbu/PyFBU.py @@ -186,17 +186,14 @@ def run(self): if self.include_gammas and nbckg > 0: gammas = [mc.Uniform('flat_gamma_{0}'.format(i), lower=0., - upper=10.) for i in range(self.nbins)] + upper=2.) for i in range(self.nbins)] gammas = mc.math.stack(gammas) # construct the Poisson constraint on gammas tau = (totalbckg/totalbckg_err)**2 - m = reciprocal(tau) - print(sqrt(m)) - print(reciprocal(sqrt(totalbckg))) gamma_poissons = [ mc.Poisson('poisson_gamma_{0}'.format(i), - mu=tau[i], observed=m[i]) for i in range(self.nbins)] + mu=gammas[i]*tau[i], observed=tau[i]) for i in range(self.nbins)] # define potential to constrain truth spectrum if self.regularization: From 961eec5bbbb43a429e2637561144589c0a4fc753 Mon Sep 17 00:00:00 2001 From: Oliver Majersky Date: Mon, 22 Jul 2019 15:40:50 +0200 Subject: [PATCH 6/7] Add ability to remove gammas for below-threshold uncert. and more configuration --- fbu/PyFBU.py | 46 ++++++++++++++++++++++++++++------------------ 1 file changed, 28 insertions(+), 18 deletions(-) diff --git a/fbu/PyFBU.py b/fbu/PyFBU.py index f2b69cd..79615c6 100644 --- a/fbu/PyFBU.py +++ b/fbu/PyFBU.py @@ -38,8 +38,9 @@ def __init__(self,data=[],response=[],background={}, self.backgroundsyst = backgroundsyst self.backgrounderr = {} self.objsyst = objsyst - self.include_gammas = False - self.gammas_treshold = 0.005 + self.include_gammas = None + self.gammas_lower = 0. + self.gammas_upper = 2. self.nbins = 0 self.systfixsigma = 0. self.smear_bckgs = {} # backgrounds to be smeared in PE (according to MC stats) @@ -64,7 +65,7 @@ def checklen(list1,list2): for bin in [self.lower,self.upper]: checklen(bin,responsetruthbins) - if self.include_gammas: + if self.include_gammas is not None: assert self.backgrounderr != {},\ 'To include gammas, must provide background MC stat uncertainties' #__________________________________________________________ @@ -103,16 +104,19 @@ def run(self): backgroundnormsysts = array([]) if nbckg>0: backgrounds = array([background[key] for key in backgroundkeys]) - # backgrounds_err_sq = array([self.backgrounderr[key]**2 for key in backgroundkeys]) - backgrounds_err_sq = array([self.backgrounderr[key] for key in backgroundkeys]) - backgrounds_err_sq = backgrounds_err_sq**2 + if self.include_gammas is not None: + backgrounds_err_sq = array([self.backgrounderr[key] for key in backgroundkeys]) + backgrounds_err_sq = backgrounds_err_sq**2 backgroundnormsysts = array([self.backgroundsyst[key] for key in backgroundkeys]) # need summed total background and it's error for gamma NPs # to take into account MC stat uncertainty of backgrounds - if self.include_gammas: + if self.include_gammas is not None: totalbckg = sum(backgrounds, axis=0) totalbckg_err = sqrt(sum(backgrounds_err_sq, axis=0)) + assert len(totalbckg) == len(self.include_gammas),\ + 'Gamma NP specification error: Inconsistent size of '\ + 'include_gammas array and the background number of bins' # unpack object systematics dictionary objsystkeys = self.objsyst['signal'].keys() @@ -184,16 +188,22 @@ def run(self): tau=1.0, **add_kwargs)) objnuisances = mc.math.stack(objnuisances) - if self.include_gammas and nbckg > 0: - gammas = [mc.Uniform('flat_gamma_{0}'.format(i), lower=0., - upper=2.) for i in range(self.nbins)] - gammas = mc.math.stack(gammas) - - # construct the Poisson constraint on gammas + if self.include_gammas is not None and nbckg > 0: + gammas = [] + gamma_poissons = [] tau = (totalbckg/totalbckg_err)**2 - gamma_poissons = [ - mc.Poisson('poisson_gamma_{0}'.format(i), - mu=gammas[i]*tau[i], observed=tau[i]) for i in range(self.nbins)] + for i,bin in enumerate(self.include_gammas): + if bin: + gammas.append(mc.Uniform('flat_gamma_{0}'.format(i), + lower=self.gammas_lower, + upper=self.gammas_upper)) + # construct the Poisson constraint on gammas + gamma_poissons.append(mc.Poisson( + 'poisson_gamma_{0}'.format(i), + mu=gammas[i]*tau[i], observed=tau[i])) + else: + gammas.append(1.) + gammas = mc.math.stack(gammas) # define potential to constrain truth spectrum if self.regularization: @@ -214,7 +224,7 @@ def unfold(): bckg = theano.dot(1. + bckgnuisances*bckgnormerr,smearedbackgrounds) - if self.include_gammas: + if self.include_gammas is not None: bckg = bckg * gammas tresmat = array(resmat) @@ -272,7 +282,7 @@ def unfold(): tmp = self.freeze_NPs[name] except KeyError as e: print('Warning: Missing NP trace', e) - if self.include_gammas: + if self.include_gammas is not None: for bin in range(self.nbins): self.nuisancestrace['gamma_{0}'.format(bin)]\ = trace['flat_gamma_{0}'.format(bin)][:] From c476bd01604fe370e4fc4cad27adb35a6180be92 Mon Sep 17 00:00:00 2001 From: Oliver Majersky Date: Fri, 6 Sep 2019 19:42:00 +0200 Subject: [PATCH 7/7] Gammas implementation with support for NP freezing for syst ranking. --- fbu/PyFBU.py | 44 +++++++++++++++++++++++++++----------------- 1 file changed, 27 insertions(+), 17 deletions(-) diff --git a/fbu/PyFBU.py b/fbu/PyFBU.py index 79615c6..64035fe 100644 --- a/fbu/PyFBU.py +++ b/fbu/PyFBU.py @@ -69,6 +69,13 @@ def checklen(list1,list2): assert self.backgrounderr != {},\ 'To include gammas, must provide background MC stat uncertainties' #__________________________________________________________ + def check_NPfrozen(self, NPname): + try: + tmp = self.freeze_NPs[NPname] + return True + except KeyError: + return False + #__________________________________________________________ def fluctuate(self, data, err=None): random.seed(self.rndseed) if err is None: @@ -193,14 +200,20 @@ def run(self): gamma_poissons = [] tau = (totalbckg/totalbckg_err)**2 for i,bin in enumerate(self.include_gammas): + NPname = 'gamma_{0}'.format(i) if bin: - gammas.append(mc.Uniform('flat_gamma_{0}'.format(i), + try: + add_kwargs['observed'] = self.freeze_NPs[NPname] + except KeyError: + add_kwargs = dict() + gammas.append(mc.Uniform('flat_' + NPname, lower=self.gammas_lower, - upper=self.gammas_upper)) + upper=self.gammas_upper, + **add_kwargs)) # construct the Poisson constraint on gammas - gamma_poissons.append(mc.Poisson( - 'poisson_gamma_{0}'.format(i), - mu=gammas[i]*tau[i], observed=tau[i])) + gamma_poissons.append(mc.Poisson('poisson_' + NPname, + mu=gammas[i]*tau[i], + observed=tau[i])) else: gammas.append(1.) gammas = mc.math.stack(gammas) @@ -278,14 +291,14 @@ def unfold(): self.nuisancestrace[name] = trace['gaus_%s'%name][:] #self.nuisancestrace[name] = copy.deepcopy(trace['gaus_%s'%name][:]) except KeyError as e: - try: - tmp = self.freeze_NPs[name] - except KeyError as e: - print('Warning: Missing NP trace', e) + if not self.check_NPfrozen(name): + print('Warning: Missing NP trace', name) if self.include_gammas is not None: for bin in range(self.nbins): - self.nuisancestrace['gamma_{0}'.format(bin)]\ - = trace['flat_gamma_{0}'.format(bin)][:] + NPname = 'gamma_{0}'.format(bin) + if self.include_gammas[bin] and not self.check_NPfrozen(NPname): + self.nuisancestrace[NPname]\ + = trace['flat_' + NPname][:] for name in objsystkeys: if self.systfixsigma==0.: try: @@ -293,12 +306,9 @@ def unfold(): self.nuisancestrace[name] = trace['flat_%s'%name][:] else: self.nuisancestrace[name] = trace['gaus_%s'%name][:] - #self.nuisancestrace[name] = copy.deepcopy(trace['gaus_%s'%name][:]) - except KeyError as e: - try: - tmp = self.freeze_NPs[name] - except KeyError as e: - print('Warning: Missing NP trace', e) + except KeyError: + if not self.check_NPfrozen(name): + print('Warning: Missing NP trace', name) if self.monitoring: from fbu import monitoring