From ac04e6a4a87489ccaf5ec8915cb9dc26b6dff3e8 Mon Sep 17 00:00:00 2001 From: chad Date: Wed, 16 Dec 2020 22:51:28 -0500 Subject: [PATCH 01/13] added n_transits to HAL function --- hawc_hal/HAL.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/hawc_hal/HAL.py b/hawc_hal/HAL.py index 569200a..faaded5 100644 --- a/hawc_hal/HAL.py +++ b/hawc_hal/HAL.py @@ -47,16 +47,19 @@ class HAL(PluginPrototype): :param flat_sky_pixels_size: size of the pixel for the flat sky projection (Hammer Aitoff) """ - def __init__(self, name, maptree, response_file, roi, flat_sky_pixels_size=0.17): + def __init__(self, name, maptree, response_file, roi, flat_sky_pixels_size=0.17, n_transits=None): # Store ROI self._roi = roi + # optionally specify n_transits + assert (n_transits==None or n_transits > 0.0), "You must specify n_transits >0" + # Set up the flat-sky projection self._flat_sky_projection = roi.get_flat_sky_projection(flat_sky_pixels_size) # Read map tree (data) - self._maptree = map_tree_factory(maptree, roi=roi) + self._maptree = map_tree_factory(maptree, roi=roi, n_transits) # Read detector response_file self._response = hawc_response_factory(response_file) From d40dfd79370489d18d24b20f91b779058253bdda Mon Sep 17 00:00:00 2001 From: chad Date: Wed, 16 Dec 2020 22:54:29 -0500 Subject: [PATCH 02/13] added transit arguments to file reading --- hawc_hal/maptree/map_tree.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/hawc_hal/maptree/map_tree.py b/hawc_hal/maptree/map_tree.py index 0b49416..b77c68b 100644 --- a/hawc_hal/maptree/map_tree.py +++ b/hawc_hal/maptree/map_tree.py @@ -17,18 +17,18 @@ import astropy.units as u -def map_tree_factory(map_tree_file, roi): +def map_tree_factory(map_tree_file, roi, n_transits=None): # Sanitize files in input (expand variables and so on) map_tree_file = sanitize_filename(map_tree_file) if os.path.splitext(map_tree_file)[-1] == '.root': - return MapTree.from_root_file(map_tree_file, roi) + return MapTree.from_root_file(map_tree_file, roi, n_transits) else: - return MapTree.from_hdf5(map_tree_file, roi) + return MapTree.from_hdf5(map_tree_file, roi, n_transits) class MapTree(object): @@ -39,9 +39,9 @@ def __init__(self, analysis_bins, roi): self._roi = roi @classmethod - def from_hdf5(cls, map_tree_file, roi): + def from_hdf5(cls, map_tree_file, roi, n_transits): - data_analysis_bins = from_hdf5_file(map_tree_file, roi) + data_analysis_bins = from_hdf5_file(map_tree_file, roi, n_transits) return cls(data_analysis_bins, roi) @@ -55,7 +55,7 @@ def from_root_file(cls, map_tree_file, roi): :return: """ - data_analysis_bins = from_root_file(map_tree_file, roi) + data_analysis_bins = from_root_file(map_tree_file, roi, n_transits) return cls(data_analysis_bins, roi) From af646954c05f45f812a78760a42c777f3edaccfd Mon Sep 17 00:00:00 2001 From: chad Date: Wed, 16 Dec 2020 22:58:27 -0500 Subject: [PATCH 03/13] added specifying transits to root analysis bins --- hawc_hal/maptree/from_root_file.py | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/hawc_hal/maptree/from_root_file.py b/hawc_hal/maptree/from_root_file.py index 63dfcfb..e84989e 100644 --- a/hawc_hal/maptree/from_root_file.py +++ b/hawc_hal/maptree/from_root_file.py @@ -41,12 +41,13 @@ def _get_bin_object(f, bin_name, suffix): return bin_tobject -def from_root_file(map_tree_file, roi): +def from_root_file(map_tree_file, roi, n_transits): """ Create a MapTree object from a ROOT file and a ROI. Do not use this directly, use map_tree_factory instead. :param map_tree_file: :param roi: + :param n_transits: :return: """ @@ -93,15 +94,17 @@ def from_root_file(map_tree_file, roi): # A transit is defined as 1 day, and totalDuration is in hours - # Get the number of transit from bin 0 (as LiFF does) - - n_transits = root_numpy.tree2array(f.Get("BinInfo"), "totalDuration") / 24.0 + # Get the number of transits from bin 0 (as LiFF does) + map_transits = root_numpy.tree2array(f.Get("BinInfo"), "totalDuration") / 24.0 # The map-maker underestimate the livetime of bins with low statistic by removing time intervals with # zero events. Therefore, the best estimate of the livetime is the maximum of n_transits, which normally # happen in the bins with high statistic - n_transits = max(n_transits) - + # Alternatively, specify n_transits + use_transits = max(map_transits) + if n_transits != None: + use_transits = n_transits + n_bins = len(data_bins_labels) # These are going to be Healpix maps, one for each data analysis bin_name @@ -145,7 +148,7 @@ def from_root_file(map_tree_file, roi): counts_hpx, bkg_hpx, active_pixels_ids=active_pixels, - n_transits=n_transits, + n_transits=use_transits, scheme='RING') else: @@ -159,7 +162,7 @@ def from_root_file(map_tree_file, roi): DenseHealpix(counts), DenseHealpix(bkg), active_pixels_ids=None, - n_transits=n_transits, + n_transits=use_transits, scheme='RING') data_analysis_bins[name] = this_data_analysis_bin From 7b09793fe64693089679b1254dc0434b49a07e20 Mon Sep 17 00:00:00 2001 From: chad Date: Wed, 16 Dec 2020 22:58:45 -0500 Subject: [PATCH 04/13] added specifying transits to hdf5 analysis bins --- hawc_hal/maptree/from_hdf5_file.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/hawc_hal/maptree/from_hdf5_file.py b/hawc_hal/maptree/from_hdf5_file.py index b894e29..609a34d 100644 --- a/hawc_hal/maptree/from_hdf5_file.py +++ b/hawc_hal/maptree/from_hdf5_file.py @@ -10,12 +10,13 @@ from .data_analysis_bin import DataAnalysisBin -def from_hdf5_file(map_tree_file, roi): +def from_hdf5_file(map_tree_file, roi, n_transits): """ Create a MapTree object from a HDF5 file and a ROI. Do not use this directly, use map_tree_factory instead. :param map_tree_file: :param roi: + :param n_transits: :return: """ @@ -88,12 +89,17 @@ def from_hdf5_file(map_tree_file, roi): # This signals the DataAnalysisBin that we are dealing with a full sky map active_pixels_user = None + # Specify n_transits (or not), default value contained in map is this_meta['n_transits'] + use_transits=this_meta['n_transits'] + if n_transits!=None: + use_transits=n_transits + # Let's now build the instance this_bin = DataAnalysisBin(bin_name, observation_hpx_map=observation_hpx_map, background_hpx_map=background_hpx_map, active_pixels_ids=active_pixels_user, - n_transits=this_meta['n_transits'], + n_transits=use_transits, scheme='RING' if this_meta['scheme'] == 0 else 'NEST') data_analysis_bins[bin_name] = this_bin From f39cb35beedc29ab76dbb67f79082b6395558ecb Mon Sep 17 00:00:00 2001 From: chad Date: Wed, 16 Dec 2020 23:35:08 -0500 Subject: [PATCH 05/13] set all arguments without keywords --- hawc_hal/HAL.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hawc_hal/HAL.py b/hawc_hal/HAL.py index faaded5..ad93587 100644 --- a/hawc_hal/HAL.py +++ b/hawc_hal/HAL.py @@ -59,7 +59,7 @@ def __init__(self, name, maptree, response_file, roi, flat_sky_pixels_size=0.17, self._flat_sky_projection = roi.get_flat_sky_projection(flat_sky_pixels_size) # Read map tree (data) - self._maptree = map_tree_factory(maptree, roi=roi, n_transits) + self._maptree = map_tree_factory(maptree, roi, n_transits) # Read detector response_file self._response = hawc_response_factory(response_file) From d62fbad32f222c3d396ed75dbca4734e85713923 Mon Sep 17 00:00:00 2001 From: chad Date: Wed, 16 Dec 2020 23:38:43 -0500 Subject: [PATCH 06/13] give n_transits to root files --- hawc_hal/maptree/map_tree.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hawc_hal/maptree/map_tree.py b/hawc_hal/maptree/map_tree.py index b77c68b..e1674be 100644 --- a/hawc_hal/maptree/map_tree.py +++ b/hawc_hal/maptree/map_tree.py @@ -46,7 +46,7 @@ def from_hdf5(cls, map_tree_file, roi, n_transits): return cls(data_analysis_bins, roi) @classmethod - def from_root_file(cls, map_tree_file, roi): + def from_root_file(cls, map_tree_file, roi, n_transits): """ Create a MapTree object from a ROOT file and a ROI. Do not use this directly, use map_tree_factory instead. From c9f42f9a45e6a37fcc29459f9223f24b67f71884 Mon Sep 17 00:00:00 2001 From: Chad Brisbois Date: Mon, 29 Mar 2021 13:46:30 -0400 Subject: [PATCH 07/13] added maptree ntransits property to allow better testing --- hawc_hal/maptree/map_tree.py | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/hawc_hal/maptree/map_tree.py b/hawc_hal/maptree/map_tree.py index e1674be..8ff6b05 100644 --- a/hawc_hal/maptree/map_tree.py +++ b/hawc_hal/maptree/map_tree.py @@ -37,11 +37,20 @@ def __init__(self, analysis_bins, roi): self._analysis_bins = analysis_bins self._roi = roi + self._n_transits = None @classmethod def from_hdf5(cls, map_tree_file, roi, n_transits): + """ + Create a MapTree object from a HDF5 file and a ROI. Do not use this directly, use map_tree_factory instead. + + :param map_tree_file: + :param roi: + :return: + """ data_analysis_bins = from_hdf5_file(map_tree_file, roi, n_transits) + self._n_transits = n_transits return cls(data_analysis_bins, roi) @@ -56,6 +65,7 @@ def from_root_file(cls, map_tree_file, roi, n_transits): """ data_analysis_bins = from_root_file(map_tree_file, roi, n_transits) + self._n_transits = n_transits return cls(data_analysis_bins, roi) @@ -96,6 +106,10 @@ def __len__(self): return len(self._analysis_bins) + @property + def n_transits(self): + return self._n_transits + @property def analysis_bins_labels(self): From fa7215c5f2dd09e43d96784ade2d3c3053dfe511 Mon Sep 17 00:00:00 2001 From: Chad Brisbois Date: Mon, 29 Mar 2021 16:38:32 -0400 Subject: [PATCH 08/13] made treatment of transits consistent betweeh hdf5 and root files, functions now return transits to pass to Maptree object --- hawc_hal/maptree/from_hdf5_file.py | 21 +++++++++++++++------ hawc_hal/maptree/from_root_file.py | 4 ++-- 2 files changed, 17 insertions(+), 8 deletions(-) diff --git a/hawc_hal/maptree/from_hdf5_file.py b/hawc_hal/maptree/from_hdf5_file.py index 92de9b4..5b0dbca 100644 --- a/hawc_hal/maptree/from_hdf5_file.py +++ b/hawc_hal/maptree/from_hdf5_file.py @@ -12,6 +12,7 @@ from ..healpix_handling import SparseHealpix, DenseHealpix from .data_analysis_bin import DataAnalysisBin +import numpy as np def from_hdf5_file(map_tree_file, roi, n_transits): """ @@ -67,6 +68,19 @@ def from_hdf5_file(map_tree_file, roi, n_transits): data_analysis_bins = collections.OrderedDict() + # Figure out the transits + transits_bins = [] + for bin_name in bin_names: + this_meta = meta_df.loc[bin_name] + # Specify n_transits (or not), default value contained in map is this_meta['n_transits'] + transits_bins.append(this_meta['n_transits']) + + # pick out the transits same as root file + use_transits = np.max(transits_bins) + if n_transits!=None: + use_transits=n_transits + + for bin_name in bin_names: this_df = analysis_bins_df.loc[bin_name] @@ -91,11 +105,6 @@ def from_hdf5_file(map_tree_file, roi, n_transits): # This signals the DataAnalysisBin that we are dealing with a full sky map active_pixels_user = None - - # Specify n_transits (or not), default value contained in map is this_meta['n_transits'] - use_transits=this_meta['n_transits'] - if n_transits!=None: - use_transits=n_transits # Let's now build the instance this_bin = DataAnalysisBin(bin_name, @@ -107,4 +116,4 @@ def from_hdf5_file(map_tree_file, roi, n_transits): data_analysis_bins[bin_name] = this_bin - return data_analysis_bins + return data_analysis_bins, use_transits diff --git a/hawc_hal/maptree/from_root_file.py b/hawc_hal/maptree/from_root_file.py index 8f48af4..58952ed 100644 --- a/hawc_hal/maptree/from_root_file.py +++ b/hawc_hal/maptree/from_root_file.py @@ -103,7 +103,7 @@ def from_root_file(map_tree_file, roi, n_transits): # zero events. Therefore, the best estimate of the livetime is the maximum of n_transits, which normally # happen in the bins with high statistic # Alternatively, specify n_transits - use_transits = max(map_transits) + use_transits = np.max(map_transits) if n_transits != None: use_transits = n_transits @@ -169,7 +169,7 @@ def from_root_file(map_tree_file, roi, n_transits): data_analysis_bins[name] = this_data_analysis_bin - return data_analysis_bins + return data_analysis_bins, use_transits def _read_partial_tree(ttree_instance, elements_to_read): From f73cde797beab93f94117d1950f08be3edf39489 Mon Sep 17 00:00:00 2001 From: Chad Brisbois Date: Mon, 29 Mar 2021 16:39:58 -0400 Subject: [PATCH 09/13] Maptree now take optional n_transits argument, if none given, pulls from data_analysis_bins --- hawc_hal/maptree/map_tree.py | 27 ++++++++++++++++++--------- 1 file changed, 18 insertions(+), 9 deletions(-) diff --git a/hawc_hal/maptree/map_tree.py b/hawc_hal/maptree/map_tree.py index 0ac3b54..f71b35a 100644 --- a/hawc_hal/maptree/map_tree.py +++ b/hawc_hal/maptree/map_tree.py @@ -35,11 +35,22 @@ def map_tree_factory(map_tree_file, roi, n_transits=None): class MapTree(object): - def __init__(self, analysis_bins, roi): + def __init__(self, analysis_bins, roi, n_transits=None): self._analysis_bins = analysis_bins self._roi = roi - self._n_transits = None + + if n_transits != None: + self._n_transits = n_transits + else: + # we need to find the transits in the map + transits=[] + # loop through bins and get the n_transits in each + for analysis_bin_key in self._analysis_bins: + transits.append(self._analysis_bins[analysis_bin_key].n_transits) + + #Get the largest + self._n_transits = np.max(np.array(transits)) @classmethod def from_hdf5(cls, map_tree_file, roi, n_transits): @@ -51,10 +62,9 @@ def from_hdf5(cls, map_tree_file, roi, n_transits): :return: """ - data_analysis_bins = from_hdf5_file(map_tree_file, roi, n_transits) - self._n_transits = n_transits - - return cls(data_analysis_bins, roi) + data_analysis_bins, transits = from_hdf5_file(map_tree_file, roi, n_transits) + + return cls(data_analysis_bins, roi, transits) @classmethod def from_root_file(cls, map_tree_file, roi, n_transits): @@ -66,10 +76,9 @@ def from_root_file(cls, map_tree_file, roi, n_transits): :return: """ - data_analysis_bins = from_root_file(map_tree_file, roi, n_transits) - self._n_transits = n_transits + data_analysis_bins, transits = from_root_file(map_tree_file, roi, n_transits) - return cls(data_analysis_bins, roi) + return cls(data_analysis_bins, roi, transits) def __iter__(self): """ From cbf54d78fb969a54266af8a424d4349e6337d8f9 Mon Sep 17 00:00:00 2001 From: Chad Brisbois Date: Mon, 29 Mar 2021 16:42:42 -0400 Subject: [PATCH 10/13] added n_transits test, checks maptree and each data_analysis_bin in it --- tests/conftest.py | 12 ++++++++++++ tests/test_n_transits.py | 34 ++++++++++++++++++++++++++++++++++ 2 files changed, 46 insertions(+) create mode 100644 tests/test_n_transits.py diff --git a/tests/conftest.py b/tests/conftest.py index f136176..e3b96d3 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -92,6 +92,18 @@ def check_map_trees(m1, m2): assert p1.nside == p2.nside assert p1.n_transits == p2.n_transits +def check_n_transits(maptree,ntransits): + # Should give maptree object and expected ntransits + + #Check maptree is assigned proper transits + assert maptree.n_transits == ntransits + + #Check each bin too + for bin_key in maptree: + + the_bin = maptree[bin_key] + assert the_bin.n_transits == ntransits + def check_responses(r1, r2): diff --git a/tests/test_n_transits.py b/tests/test_n_transits.py new file mode 100644 index 0000000..be9a130 --- /dev/null +++ b/tests/test_n_transits.py @@ -0,0 +1,34 @@ +from hawc_hal.maptree.map_tree import map_tree_factory +from hawc_hal import HealpixConeROI +import os +import pytest +from conftest import check_n_transits + + +def test_transits(maptree, + roi): + + roi_ = roi + + #Case 1: not specifying transits + + # Dont specify n_transits, use value in maptree + maptree_none = map_tree_factory(maptree, roi_) + + n_transits = maptree_none.n_transits + + #This is somewhat trivial, but something has definitely gone wrong if this fails + check_n_transits(maptree_none,n_transits) + + + #Case 2: specify transits + + # specify transits, in this case 777.7 + n_transits_specify = 777.7 + + # specify n_transits + maptree_specify = map_tree_factory(maptree, roi_, n_transits_specify) + + # Does maptree return specified transits? + check_n_transits(maptree_specify,n_transits_specify) + From 60d4eebfd7300cba921087c9a5841f9d120ce976 Mon Sep 17 00:00:00 2001 From: Chad Brisbois Date: Mon, 29 Mar 2021 17:05:58 -0400 Subject: [PATCH 11/13] unharshed the error message, no longer crashes and gives user info about behavior --- hawc_hal/HAL.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/hawc_hal/HAL.py b/hawc_hal/HAL.py index d25359f..44b4fea 100644 --- a/hawc_hal/HAL.py +++ b/hawc_hal/HAL.py @@ -50,13 +50,19 @@ class HAL(PluginPrototype): :param flat_sky_pixels_size: size of the pixel for the flat sky projection (Hammer Aitoff) """ - def __init__(self, name, maptree, response_file, roi, flat_sky_pixels_size=0.17, n_transits=None): + def __init__(self, name, maptree, response_file, roi, flat_sky_pixels_size=0.17, set_transits=-1): # Store ROI self._roi = roi # optionally specify n_transits - assert (n_transits==None or n_transits > 0.0), "You must specify n_transits >0" + if (num_transits > 0.0): + log.info("Setting Transits to {0}".format(set_transits)) + n_transits=set_transits + else: + log.info("Using transits contained in the maptree") + n_transits=None + # Set up the flat-sky projection self._flat_sky_projection = roi.get_flat_sky_projection(flat_sky_pixels_size) From 297fe998cd598e62d037e0afa17a96ed255974de Mon Sep 17 00:00:00 2001 From: Chad Brisbois Date: Mon, 29 Mar 2021 22:21:26 -0400 Subject: [PATCH 12/13] fix Linter error --- hawc_hal/HAL.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hawc_hal/HAL.py b/hawc_hal/HAL.py index 44b4fea..5a56609 100644 --- a/hawc_hal/HAL.py +++ b/hawc_hal/HAL.py @@ -56,7 +56,7 @@ def __init__(self, name, maptree, response_file, roi, flat_sky_pixels_size=0.17, self._roi = roi # optionally specify n_transits - if (num_transits > 0.0): + if (set_transits > 0.0): log.info("Setting Transits to {0}".format(set_transits)) n_transits=set_transits else: From 0e3238b6a3a3f59dc873ad69afdf4a8d2d5cc708 Mon Sep 17 00:00:00 2001 From: Chad Brisbois Date: Mon, 29 Mar 2021 22:34:50 -0400 Subject: [PATCH 13/13] access transits array directly in hdf5 files, symmetrize if transit statements for root and hdf5 --- hawc_hal/maptree/from_hdf5_file.py | 16 +++++----------- hawc_hal/maptree/from_root_file.py | 12 +++++++----- 2 files changed, 12 insertions(+), 16 deletions(-) diff --git a/hawc_hal/maptree/from_hdf5_file.py b/hawc_hal/maptree/from_hdf5_file.py index 5b0dbca..925cd7b 100644 --- a/hawc_hal/maptree/from_hdf5_file.py +++ b/hawc_hal/maptree/from_hdf5_file.py @@ -64,21 +64,15 @@ def from_hdf5_file(map_tree_file, roi, n_transits): bin_names = analysis_bins_df.index.levels[0] - # Loop over them and build a DataAnalysisBin instance for each one - data_analysis_bins = collections.OrderedDict() # Figure out the transits - transits_bins = [] - for bin_name in bin_names: - this_meta = meta_df.loc[bin_name] - # Specify n_transits (or not), default value contained in map is this_meta['n_transits'] - transits_bins.append(this_meta['n_transits']) - # pick out the transits same as root file - use_transits = np.max(transits_bins) - if n_transits!=None: - use_transits=n_transits + if n_transits is None: + use_transits = meta_df["n_transits"].max() + else: + # specified n_transits + use_transits = n_transits for bin_name in bin_names: diff --git a/hawc_hal/maptree/from_root_file.py b/hawc_hal/maptree/from_root_file.py index 58952ed..ab8f0cd 100644 --- a/hawc_hal/maptree/from_root_file.py +++ b/hawc_hal/maptree/from_root_file.py @@ -95,16 +95,18 @@ def from_root_file(map_tree_file, roi, n_transits): data_bins_labels = [ str(i) for i in data_bins_labels ] - # A transit is defined as 1 day, and totalDuration is in hours - # Get the number of transits from bin 0 (as LiFF does) - map_transits = root_numpy.tree2array(f.Get("BinInfo"), "totalDuration") / 24.0 # The map-maker underestimate the livetime of bins with low statistic by removing time intervals with # zero events. Therefore, the best estimate of the livetime is the maximum of n_transits, which normally # happen in the bins with high statistic # Alternatively, specify n_transits - use_transits = np.max(map_transits) - if n_transits != None: + + if n_transits is None: + # A transit is defined as 1 day, and totalDuration is in hours + # Get the number of transits from bin 0 (as LiFF does) + map_transits = root_numpy.tree2array(f.Get("BinInfo"), "totalDuration") / 24.0 + use_transits = np.max(map_transits) + else: use_transits = n_transits n_bins = len(data_bins_labels)