Skip to content

Commit 04aee41

Browse files
authored
Implementation of diffuse flux in the background (#449)
* make SoB splines class * 1. add bkg spatial into SoB_splines class 2. remove kwargs from class init, instead takes a dict where name, smoothing order, and gamma precision should be in. NO default vals! The latter 2 can be of str type. 3. bkg_weights class method returns NotImplementedError instead of array of 1s 4. remove obsolete environment smoothing & precision keys 5. remove redundant if clauses pertaining to making smoothing order an int & environment key in make_2d_spline_from_hist method 6. type annotate smoothing order to int where needed 7. remove kwargs from methods, use the class params instead 8. change the plot paths in make_individual_spline_set method to include a suffix with the subclass name if not "no_difffuse". The sig, SoB hist & splines plots will now be under the <gamma>/precision<float>_smoothing<int> dir 9. change the SoB_spline_path func args to include the sob dict as kwargs in make_splines & load_spline methods 10. BugFix: use per-flavour fluxes for SPL & BPL subclasses 11. BugFix: make the diffuse + atm weights an np.array 12. BPL subclass get_diffuse_flux method returns fluxes and respective masks for the 2 energy ranges. That accommodates the subclass bkg_weights method to comply with returning a single weights array _Note to self: commit more oft!_ * import SoB_splines class to use its staticmethod make_plot * add SoB class name in the spline paths if it's NOT the default "no_difffuse" (ie atm-only). Adds extra arg in the smoothing_precision_string, llh_energy_hash_pickles, and bkg_spline_path func * add simulate_bkg_with_diffuse method in the NTSeason. Same as the simulate_background(), but takes a bkg_weights func as an arg, so that the diffuse + atm weights are calculated for given MC * 1. add simulate_bkg_with_diffuse() in Season, which raises NotImplementedError (only for NT!!) 2. add SoB arg in make_background_spatial(). The SoB_splines.mmake_background_spline() is then called 3. change imports * Add an optional SoB arg when creating SpatialPDF. This is only relevant for the BackgroundSpatialPDF, particularly for the ZenithSpline where SoB_splines class is needed for the create_background_function method. For UniformSolidAngle bkg spatial pdf the SoB = None * Changes in BaseInjector: 1. when initializing, make sob dict with default vals for bkg_model_name, smoothing_order, and gamma_precision and instantiate SoB_splines. 2. pass the SoB class when creating the spatial_pdf if the "bkg_spatial_pdf" DOES NOT exist in the injection_spatial_pdf dict (BackgroundSpatialPDF defaults to ZenithSpline), otherwise is None NOTE: not optimal and subject to change, since it depends on BackgroundSpatialPDF.create method!!! 3. when creating the dataset, call season.simulate_bkg_with_diffuse() w/ Sob_splines.bkg_weights as its Callable arg if sob_name isn't "no_difffuse", otherwise use the regular simulate_background() * BUGFIX: when creating dataset in BaseInjector and SoB spline name isn't "no_difffuse", have a try except clause to catch cases when datasets other than NT are used. If that's the case a NotImplementedError is raised and simulate_background() is used instead. Note: currently simulate_bkg_with_diffuse() is only implemented for NTSeason but can be extended to all datasets w/ MC * changes to LLH 1. when initializing LLH a) make sob dict with default vals for smoothing, precision, and SoB name. If sob dict isn't passed in the llh dict, then it's created with smoothing and precision vals taken from the llh dict for backwards compatibility. Create SoB spline and pass its smoothing_order and gamma_precision properties to the respective LLH ones. b) when creating the SpatialPDF, pass SoB instance if bkg spatial is not specified in the llh spatial dict, otherwise pass None. !!Hacky method, dependent on the hardcoded default "zenith spline" for the bkg spatial. Prone to error if that changes!!!! 2. BUGFIX: update return val in SpatialLLH create_energy_function. Calls self.spatial_pdf.background_spatial instead of obsolete property 3. changes in FixedEnergyLLH: a) delete redundant smoothing order and precision properties assignment when initializing b) update llh_energy_hash_pickles() c) replace old methods to SoB class methods where needed d) BUGFIX: dump proper sinDEC bins from season when making SoB pickle in create_energy_weighting_function() 4. changes in StandardLLH a) replace old methods to SoB class methods where needed b) update SoB_spline_path() 5. BUGFIX: in generate_dynamic_flare_class() raise NotImplementedError if mh_name is either "spatial", or "fixed_energy" since there is no SoB energy cache in these 2 * add warning logging when initializing minimizer if the SoB dicts in the llh & inj dicts are different * delete redundant warning log in the minimizer for when running fixed gamma trials with different gamma in the llh dict than the one in inj dict * BUGFIX: have a try except clause in the submit method of Submitter for when LowMemoryInjector is used and the make_band_mask needs to be called. Fixes cases when you submit jobs on the cluster that don't include inj dict * BUGFIX: change MockUnblindedInjector to return scrambled data for all seasons. Substitute load_background_model() with use_data_for_trials() so the exp data are correctly loaded when initializing. Bypass the season's simulate_background() in create_dataset() by basically copying the Season.simulate_background() method before applying the angular_error_modifier (bad practice but can't think of sth better) * update public all_sky_3_year module used in tests. Make default SoB class (no_difffuse, flarestack smoothing order and precision) to call make_individual_spline_set method * fix linter: BPLDiffuseSpline inherits from base class and not SPLDiffuseSpline; return correct types for masks in its get_diffuse_flux method * isort fix * exit ResultsHandler after making TS and params plots if only bkg trials
1 parent 435a6a4 commit 04aee41

12 files changed

Lines changed: 880 additions & 566 deletions

File tree

flarestack/cluster/submitter.py

Lines changed: 6 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -95,8 +95,12 @@ def submit(self, mh_dict):
9595
if self.remove_old_results:
9696
self._clean_injection_values_and_pickled_results(self.mh_dict["name"])
9797
if self.use_cluster:
98-
if mh_dict["inj_dict"]["injector_name"] == "low_memory_injector":
99-
make_band_mask(mh_dict=copy.deepcopy(mh_dict))
98+
try:
99+
inj_name = mh_dict["inj_dict"]["injector_name"]
100+
if inj_name == "low_memory_injector":
101+
make_band_mask(mh_dict=copy.deepcopy(mh_dict))
102+
except KeyError:
103+
pass
100104

101105
self.submit_cluster(mh_dict)
102106
else:

flarestack/core/injector.py

Lines changed: 47 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -15,6 +15,7 @@
1515
from flarestack.core.time_pdf import TimePDF, read_t_pdf_dict
1616
from flarestack.shared import band_mask_cache_name, k_to_flux
1717
from flarestack.utils.catalogue_loader import calculate_source_weight
18+
from flarestack.utils.make_SoB_splines import SoB_splines
1819

1920
if TYPE_CHECKING:
2021
from flarestack.data import Season, SeasonWithMC
@@ -87,14 +88,35 @@ def __init__(self, season: "Season", sources: Table, **kwargs) -> None:
8788
if len(sources) > 0:
8889
self.weight_scale = calculate_source_weight(self.sources)
8990

91+
# get SoB splines dict
92+
try:
93+
sob_dict = kwargs["sob_dict"]
94+
except KeyError:
95+
sob_dict = {}
96+
97+
sob_dict.setdefault("smoothing_order", "flarestack")
98+
sob_dict.setdefault("gamma_precision", "flarestack")
99+
self.sob_name = sob_dict.setdefault("bkg_model_name", "no_difffuse")
100+
self.sob = SoB_splines.create(season, sob_dict)
101+
90102
try:
91103
self.sig_time_pdf = TimePDF.create(
92104
kwargs["injection_sig_time_pdf"], season.get_time_pdf()
93105
)
94106
# self.bkg_time_pdf = TimePDF.create(kwargs["injection_bkg_time_pdf"],
95107
# season.get_time_pdf())
96108
self.energy_pdf = EnergyPDF.create(kwargs["injection_energy_pdf"])
97-
self.spatial_pdf = SpatialPDF(kwargs["injection_spatial_pdf"], season)
109+
110+
# bkg_spatial_pdf = zenith_spline if not specified
111+
# for this case SoB dict is needed
112+
if "bkg_spatial_pdf" not in kwargs["injection_spatial_pdf"].keys():
113+
self.spatial_pdf = SpatialPDF(
114+
kwargs["injection_spatial_pdf"], season, self.sob
115+
)
116+
else:
117+
self.spatial_pdf = SpatialPDF(
118+
kwargs["injection_spatial_pdf"], season, None
119+
)
98120
except KeyError:
99121
raise Exception(
100122
"Injection Arguments missing. \n "
@@ -166,9 +188,19 @@ def create_dataset(
166188
:param angular_error_modifier: AngularErrorModifier to change angular errors
167189
:return: Simulated dataset
168190
"""
169-
bkg_events, n_excluded = self.season.simulate_background(
170-
self.sources, self.spatial_box_width
171-
)
191+
if self.sob_name != "no_difffuse":
192+
try:
193+
bkg_events, n_excluded = self.season.simulate_bkg_with_diffuse(
194+
self.sources, self.spatial_box_width, self.sob.bkg_weights
195+
)
196+
except NotImplementedError:
197+
bkg_events, n_excluded = self.season.simulate_background(
198+
self.sources, self.spatial_box_width
199+
)
200+
else:
201+
bkg_events, n_excluded = self.season.simulate_background(
202+
self.sources, self.spatial_box_width
203+
)
172204

173205
if scale > 0.0:
174206
sig_events = self.inject_signal(scale)
@@ -729,8 +761,7 @@ class MockUnblindedInjector:
729761

730762
def __init__(self, season: "Season", sources=np.nan, **kwargs):
731763
self.season = season
732-
self._raw_data = season.get_exp_data()
733-
season.load_background_model()
764+
self.season.use_data_for_trials()
734765

735766
def create_dataset(
736767
self, scale: float, angular_error_modifier=None
@@ -742,11 +773,18 @@ def create_dataset(
742773
seed = int(123456)
743774
np.random.seed(seed)
744775

745-
simulated_data, n_excluded = self.season.simulate_background(Table(), None)
776+
# copy Season.simulate_background()
777+
scrambled_data = self.season.pseudo_background()
778+
if self.season._subselection_fraction is not None:
779+
scrambled_data = np.random.choice(
780+
scrambled_data,
781+
int(len(scrambled_data) * self.season._subselection_fraction),
782+
)
783+
746784
if angular_error_modifier is not None:
747-
simulated_data = angular_error_modifier.pull_correct_static(simulated_data)
785+
scrambled_data = angular_error_modifier.pull_correct_static(scrambled_data)
748786

749-
return simulated_data, n_excluded
787+
return scrambled_data, 0
750788

751789

752790
class TrueUnblindedInjector:

flarestack/core/llh.py

Lines changed: 38 additions & 49 deletions
Original file line numberDiff line numberDiff line change
@@ -20,20 +20,13 @@
2020
from flarestack.shared import (
2121
SoB_spline_path,
2222
acceptance_path,
23-
default_gamma_precision,
24-
default_smoothing_order,
2523
llh_energy_hash_pickles,
2624
)
2725
from flarestack.utils.create_acceptance_functions import (
2826
dec_range,
2927
make_acceptance_season,
3028
)
31-
from flarestack.utils.make_SoB_splines import (
32-
create_2d_ratio_hist,
33-
load_spline,
34-
make_2d_spline_from_hist,
35-
make_individual_spline_set,
36-
)
29+
from flarestack.utils.make_SoB_splines import SoB_splines, get_gamma_precision
3730

3831
logger = logging.getLogger(__name__)
3932

@@ -100,7 +93,26 @@ def __init__(self, season, sources, llh_dict):
10093
self.season = season
10194
self.sources = sources
10295
self.llh_dict = llh_dict
103-
self.spatial_pdf = SpatialPDF(llh_dict["llh_spatial_pdf"], season)
96+
97+
try:
98+
sob_dict = llh_dict["sob_dict"]
99+
sob_dict.setdefault("smoothing_order", "flarestack")
100+
sob_dict.setdefault("gamma_precision", "flarestack")
101+
except KeyError:
102+
sob_dict = dict()
103+
smoothing = llh_dict.get("smoothing_order", "flarestack")
104+
precision = llh_dict.get("gamma_precision", "flarestack")
105+
sob_dict = {"smoothing_order": smoothing, "gamma_precision": precision}
106+
107+
self.sob_name = sob_dict.setdefault("bkg_model_name", "no_difffuse")
108+
self.sob = SoB_splines.create(season, sob_dict)
109+
self.smoothing_order = self.sob.smoothing_order
110+
self.precision = get_gamma_precision(self.sob.gamma_precision)
111+
112+
if "bkg_spatial_pdf" not in llh_dict["llh_spatial_pdf"].keys():
113+
self.spatial_pdf = SpatialPDF(llh_dict["llh_spatial_pdf"], season, self.sob)
114+
else:
115+
self.spatial_pdf = SpatialPDF(llh_dict["llh_spatial_pdf"], season, None)
104116
self.spatial_box_width = llh_dict.get(
105117
"spatial_box_width", default_spacial_box_width
106118
)
@@ -423,7 +435,7 @@ def create_energy_function(self):
423435
del exp
424436

425437
# return lambda x: data_rate
426-
return lambda x: np.exp(self.bkg_spatial(np.sin(x))) * data_rate
438+
return lambda x: self.spatial_pdf.background_spatial(x) * data_rate
427439

428440
def create_llh_function(
429441
self, data: Table, n_excluded: int, pull_corrector, weight_f=None
@@ -522,22 +534,6 @@ def __init__(self, season, sources, llh_dict):
522534
"again."
523535
)
524536

525-
# defines the order of the splines used in the building of the energy PDF
526-
self.smoothing_order = None
527-
smoothing_order = llh_dict.get("smoothing_order", "flarestack")
528-
if isinstance(smoothing_order, str):
529-
self.smoothing_order = default_smoothing_order[smoothing_order]
530-
else:
531-
self.smoothing_order = smoothing_order
532-
533-
# used to construct the support points in gamma when the energy PDF is built
534-
self.precision = None
535-
precision = llh_dict.get("gamma_precision", "flarestack")
536-
if isinstance(precision, str):
537-
self.precision = default_gamma_precision[precision]
538-
else:
539-
self.precision = precision
540-
541537
LLH.__init__(self, season, sources, llh_dict)
542538

543539
def create_energy_functions(self):
@@ -548,7 +544,9 @@ def create_energy_functions(self):
548544
:return: Acceptance function, energy_weighting_function
549545
"""
550546

551-
SoB_path, acc_path = llh_energy_hash_pickles(self.llh_dict, self.season)
547+
SoB_path, acc_path = llh_energy_hash_pickles(
548+
self.llh_dict, self.season, self.sob_name
549+
)
552550

553551
# Set up acceptance function, creating values if they have not been
554552
# created before
@@ -577,7 +575,7 @@ def acc_f(source, params=None):
577575
with open(SoB_path, "rb") as f:
578576
[dec_vals, ratio_hist] = pickle.load(f)
579577

580-
spline = make_2d_spline_from_hist(
578+
spline = self.sob.make_2d_spline_from_hist(
581579
np.array(ratio_hist), dec_vals, self.season.log_e_bins, self.smoothing_order
582580
)
583581

@@ -636,12 +634,12 @@ def create_energy_weighting_function(self, SoB_path):
636634

637635
# dec_range = self.season["sinDec bins"]
638636

639-
ratio_hist = create_2d_ratio_hist(
637+
ratio_hist = self.sob.create_2d_ratio_hist(
640638
exp=self.season.get_background_model(),
641639
mc=self.season.get_pseudo_mc(),
642-
sin_dec_bins=dec_range,
640+
sin_dec_bins=self.season.sin_dec_bins,
643641
log_e_bins=self.season.log_e_bins,
644-
weight_function=self.energy_pdf.weight_mc,
642+
weight_f=self.energy_pdf.weight_mc,
645643
)
646644

647645
try:
@@ -652,7 +650,7 @@ def create_energy_weighting_function(self, SoB_path):
652650
logger.info("Saving to {0}".format(SoB_path))
653651

654652
with open(SoB_path, "wb") as f:
655-
pickle.dump([dec_range, ratio_hist], f)
653+
pickle.dump([self.season.sin_dec_bins, ratio_hist], f)
656654

657655
def create_kwargs(
658656
self,
@@ -742,11 +740,7 @@ def __init__(self, season, sources, llh_dict):
742740
# Sets precision for energy SoB
743741
# self.precision = .1
744742

745-
self.SoB_spline_2Ds = load_spline(
746-
self.season,
747-
smoothing_order=self.smoothing_order,
748-
gamma_precision=self.precision,
749-
)
743+
self.SoB_spline_2Ds = self.sob.load_spline()
750744

751745
if self.SoB_spline_2Ds:
752746
logger.debug("Loaded {0} splines.".format(len(self.SoB_spline_2Ds)))
@@ -778,11 +772,7 @@ def create_energy_functions(self):
778772
:return: Acceptance function, energy_weighting_function
779773
"""
780774

781-
SoB_path = SoB_spline_path(
782-
self.season,
783-
smoothing_order=self.smoothing_order,
784-
gamma_precision=self.precision,
785-
)
775+
SoB_path = SoB_spline_path(self.season, **self.sob.sob_dict)
786776
acc_path = acceptance_path(self.season)
787777

788778
# Set up acceptance function, creating values if they have not been
@@ -807,12 +797,7 @@ def create_energy_functions(self):
807797
# Checks if energy weighting functions have been created
808798

809799
if not os.path.isfile(SoB_path):
810-
make_individual_spline_set(
811-
self.season,
812-
SoB_path,
813-
smoothing_order=self.smoothing_order,
814-
gamma_precision=self.precision,
815-
)
800+
self.sob.make_individual_spline_set(self.season, SoB_path)
816801

817802
return acc_f, None
818803

@@ -1708,8 +1693,12 @@ def generate_dynamic_flare_class(season, sources, llh_dict):
17081693
except KeyError:
17091694
raise KeyError("No LLH specified.")
17101695

1711-
# Set up dynamic inheritance
1696+
if mh_name in ["spatial", "fixed_energy"]:
1697+
raise NotImplementedError(
1698+
"Need an LLh that creates SoB energy cache, choose another"
1699+
)
17121700

1701+
# Set up dynamic inheritance
17131702
try:
17141703
ParentLLH = LLH.subclasses[mh_name]
17151704
except KeyError:

flarestack/core/minimisation.py

Lines changed: 12 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -222,10 +222,18 @@ def __init__(self, mh_dict):
222222
"gamma" in self.llh_dict["llh_energy_pdf"].keys()
223223
), "Running trials with fixed gamma but no gamma passed in the llh energy pdf"
224224
self.llh_gamma = self.llh_dict["llh_energy_pdf"]["gamma"]
225-
if self.llh_gamma != self.inj_dict["injection_energy_pdf"]["gamma"]:
226-
logger.warning(
227-
f"Fixing gamma to {self.llh_gamma} for llh but injection is with {self.inj_dict['injection_energy_pdf']['gamma']}"
228-
)
225+
226+
# check if SoB spline dict is same in the llh & inj dicts
227+
if (
228+
"sob_dict" in self.inj_dict
229+
and "sob_dict" in self.llh_dict
230+
and self.llh_dict["sob_dict"] != self.inj_dict["sob_dict"]
231+
):
232+
logger.warning(
233+
"Passed 'sob_dict' aren't the same for injector and llh: "
234+
+ f"inj = {self.inj_dict['sob_dict']}, llh = {self.llh_dict['sob_dict']}"
235+
+ " Proceed with caution!"
236+
)
229237

230238
p0, bounds, names = self.return_parameter_info(mh_dict)
231239

flarestack/core/results.py

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -141,6 +141,7 @@ def __init__(
141141
# just make the TS and param distr plots
142142
if len(self.scale_values) == 1 and self.scale_values[0] == 0:
143143
self.make_plots(self.scale_labels[0])
144+
return
144145

145146
# Create fit bias plots
146147
# this expects flux_to_ns to be set

0 commit comments

Comments
 (0)