From c155fa8c1b6875cfa9c45e60553c9ec938f74c8f Mon Sep 17 00:00:00 2001 From: YL <20972148+LiYunyang@users.noreply.github.com> Date: Sun, 1 Dec 2024 22:30:38 -0600 Subject: [PATCH 1/4] [perf] reproj partial maps - option to skip map pixels in ``ReprojectMaps`` --- .gitignore | 1 + maps/python/map_modules.py | 48 ++++++++++++++++++++++++++++++++++---- maps/src/maputils.cxx | 16 ++++++++++--- 3 files changed, 58 insertions(+), 7 deletions(-) diff --git a/.gitignore b/.gitignore index 808b3d1a..25402965 100644 --- a/.gitignore +++ b/.gitignore @@ -13,3 +13,4 @@ doc/module-autodocs.rst doc/moddoc_* *.g3 .DS_Store +.idea/ \ No newline at end of file diff --git a/maps/python/map_modules.py b/maps/python/map_modules.py index cd7b0a31..dc456030 100644 --- a/maps/python/map_modules.py +++ b/maps/python/map_modules.py @@ -871,15 +871,26 @@ class ReprojectMaps(object): weighted : bool If True (default), ensure that maps have had weights applied before reprojection. Otherwise, reproject maps without checking the weights. + partial : bool=False + If True, the reproj will be performed on a partial map (of the output map), + defined by the mask. If the mask is not provided, it will be detrmined from + the non-zero pixels of the first reprojected map. + mask : G3SkyMapMask, G3SkyMap, or np.ndarray, Optional. + Mask to be used for partial reproject. This should be of the same size as the output map. + For numpy array, all zeros/inf/nan/hp.UNSEEN pixels are skipped. """ - def __init__(self, map_stub=None, rebin=1, interp=False, weighted=True): + def __init__(self, map_stub=None, rebin=1, interp=False, weighted=True, partial=False, mask=None): assert map_stub is not None, "map_stub argument required" self.stub = map_stub.clone(False) self.stub.pol_type = None self.rebin = rebin self.interp = interp self.weighted = weighted + self._mask = None + self.partial = partial + + self.mask = mask def __call__(self, frame): if isinstance(frame, core.G3Frame) and frame.type != core.G3FrameType.Map: @@ -905,15 +916,44 @@ def __call__(self, frame): if key in "TQUH": mnew = self.stub.clone(False) - maps.reproj_map(m, mnew, rebin=self.rebin, interp=self.interp) + maps.reproj_map(m, mnew, rebin=self.rebin, interp=self.interp, mask=self.mask) elif key in ["Wpol", "Wunpol"]: mnew = maps.G3SkyMapWeights(self.stub) for wkey in mnew.keys(): maps.reproj_map( - m[wkey], mnew[wkey], rebin=self.rebin, interp=self.interp + m[wkey], mnew[wkey], rebin=self.rebin, interp=self.interp, mask=self.mask ) frame[key] = mnew - + self.mask = mnew return frame + + @property + def mask(self): + return self._mask + + @mask.setter + def mask(self, mask): + if mask is None: + return + if self._mask is None and self.partial: + if isinstance(mask, maps.G3SkyMapMask): + self._mask = mask + elif isinstance(mask, maps.G3SkyMap): + self._mask = maps.G3SkyMapMask(mask, use_data=True, zero_nans=True, zero_infs=True) + elif isinstance(mask, np.ndarray): + from healpy import UNSEEN + tmp = self.stub.clone(False) + mask_copy = np.ones(mask.shape, dtype=int) + bad = np.logical_or.reduce([ + np.isnan(mask), + np.isinf(mask), + mask==0, + mask==UNSEEN + ]) + mask_copy[bad] = 0 + tmp[:] = mask_copy + self._mask = maps.G3SkyMapMask(tmp, use_data=True) + else: + raise TypeError("Mask must be a G3SkyMapMask, G3SkyMap, or numpy array") diff --git a/maps/src/maputils.cxx b/maps/src/maputils.cxx index eddabe72..e8b12df3 100644 --- a/maps/src/maputils.cxx +++ b/maps/src/maputils.cxx @@ -279,7 +279,8 @@ void FlattenPol(FlatSkyMapPtr Q, FlatSkyMapPtr U, G3SkyMapWeightsPtr W, double h } -void ReprojMap(G3SkyMapConstPtr in_map, G3SkyMapPtr out_map, int rebin, bool interp) +void ReprojMap(G3SkyMapConstPtr in_map, G3SkyMapPtr out_map, int rebin, bool interp, + G3SkyMapMaskConstPtr out_map_mask) { bool rotate = false; // no transform Quat q_rot; // quaternion for rotating from output to input coordinate system @@ -310,8 +311,13 @@ void ReprojMap(G3SkyMapConstPtr in_map, G3SkyMapPtr out_map, int rebin, bool int out_map->pol_conv = in_map->pol_conv; } + size_t stop = out_map->size(); if (rebin > 1) { - for (size_t i = 0; i < out_map->size(); i++) { + for (size_t i = 0; i < stop; i++) { + if (!!out_map_mask && !out_map_mask->at(i)) { + (*out_map)[i] = 0; + continue; + } double val = 0; auto quats = out_map->GetRebinQuats(i, rebin); if (rotate) @@ -328,7 +334,11 @@ void ReprojMap(G3SkyMapConstPtr in_map, G3SkyMapPtr out_map, int rebin, bool int } } } else { - for (size_t i = 0; i < out_map->size(); i++) { + for (size_t i = 0; i < stop; i++) { + if (!!out_map_mask && !out_map_mask->at(i)) { + (*out_map)[i] = 0; + continue; + } double val = 0; auto q = out_map->PixelToQuat(i); if (rotate) From 42537ffa3a1f17400c0d4f2a519b79854582ede2 Mon Sep 17 00:00:00 2001 From: YL <20972148+LiYunyang@users.noreply.github.com> Date: Mon, 2 Dec 2024 10:52:54 -0600 Subject: [PATCH 2/4] [doc] typo --- maps/python/map_modules.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/maps/python/map_modules.py b/maps/python/map_modules.py index dc456030..e9d5c1b4 100644 --- a/maps/python/map_modules.py +++ b/maps/python/map_modules.py @@ -873,7 +873,7 @@ class ReprojectMaps(object): reprojection. Otherwise, reproject maps without checking the weights. partial : bool=False If True, the reproj will be performed on a partial map (of the output map), - defined by the mask. If the mask is not provided, it will be detrmined from + defined by the mask. If the mask is not provided, it will be determined from the non-zero pixels of the first reprojected map. mask : G3SkyMapMask, G3SkyMap, or np.ndarray, Optional. Mask to be used for partial reproject. This should be of the same size as the output map. From 3c3085068127518919f4eab0abbfa43a1ffc9e6e Mon Sep 17 00:00:00 2001 From: YL <20972148+LiYunyang@users.noreply.github.com> Date: Mon, 2 Dec 2024 16:09:13 -0600 Subject: [PATCH 3/4] [fix] missing py binding --- maps/python/map_modules.py | 15 ++++++++++----- maps/src/maputils.cxx | 6 ++++-- 2 files changed, 14 insertions(+), 7 deletions(-) diff --git a/maps/python/map_modules.py b/maps/python/map_modules.py index e9d5c1b4..872ae7d3 100644 --- a/maps/python/map_modules.py +++ b/maps/python/map_modules.py @@ -880,7 +880,8 @@ class ReprojectMaps(object): For numpy array, all zeros/inf/nan/hp.UNSEEN pixels are skipped. """ - def __init__(self, map_stub=None, rebin=1, interp=False, weighted=True, partial=False, mask=None): + def __init__(self, map_stub=None, rebin=1, interp=False, weighted=True, + partial=False, mask=None): assert map_stub is not None, "map_stub argument required" self.stub = map_stub.clone(False) self.stub.pol_type = None @@ -916,13 +917,15 @@ def __call__(self, frame): if key in "TQUH": mnew = self.stub.clone(False) - maps.reproj_map(m, mnew, rebin=self.rebin, interp=self.interp, mask=self.mask) + maps.reproj_map(m, mnew, rebin=self.rebin, interp=self.interp, + mask=self.mask) elif key in ["Wpol", "Wunpol"]: mnew = maps.G3SkyMapWeights(self.stub) for wkey in mnew.keys(): maps.reproj_map( - m[wkey], mnew[wkey], rebin=self.rebin, interp=self.interp, mask=self.mask + m[wkey], mnew[wkey], rebin=self.rebin, interp=self.interp, + mask=self.mask ) frame[key] = mnew @@ -941,7 +944,8 @@ def mask(self, mask): if isinstance(mask, maps.G3SkyMapMask): self._mask = mask elif isinstance(mask, maps.G3SkyMap): - self._mask = maps.G3SkyMapMask(mask, use_data=True, zero_nans=True, zero_infs=True) + self._mask = maps.G3SkyMapMask(mask, use_data=True, zero_nans=True, + zero_infs=True) elif isinstance(mask, np.ndarray): from healpy import UNSEEN tmp = self.stub.clone(False) @@ -956,4 +960,5 @@ def mask(self, mask): tmp[:] = mask_copy self._mask = maps.G3SkyMapMask(tmp, use_data=True) else: - raise TypeError("Mask must be a G3SkyMapMask, G3SkyMap, or numpy array") + raise TypeError("Mask must be a G3SkyMapMask, G3SkyMap, " + "or numpy array") diff --git a/maps/src/maputils.cxx b/maps/src/maputils.cxx index e8b12df3..b7cee8d0 100644 --- a/maps/src/maputils.cxx +++ b/maps/src/maputils.cxx @@ -582,7 +582,8 @@ PYBINDINGS("maps") "the appropriate rotation to the Q and u elements of the associated weights."); bp::def("reproj_map", ReprojMap, - (bp::arg("in_map"), bp::arg("out_map"), bp::arg("rebin")=1, bp::arg("interp")=false), + (bp::arg("in_map"), bp::arg("out_map"), bp::arg("rebin")=1, bp::arg("interp")=false, + bp::arg("mask")=bp::object()), "Reprojects the data from in_map onto out_map. out_map can have a different " "projection, size, resolution, etc. Optionally account for sub-pixel " "structure by setting rebin > 1 and/or enable bilinear interpolation of " @@ -590,7 +591,8 @@ PYBINDINGS("maps") "attributes to rotate between Equatorial and Galactic coordinate systems. " "Use the maps' pol_conv attributes to switch between COSMO and IAU " "polarization conventions. If output attributes are not set, they will be " - "copied from the input map."); + "copied from the input map. out_map_mask, if given, skip the unused pixels" + "and set these pixels to 0."); bp::def("get_map_moments", GetMapMoments, (bp::arg("map"), bp::arg("mask")=G3SkyMapMaskConstPtr(), bp::arg("order")=2, From b1b117ef8e6899945969346b7fe80d53732b6ec8 Mon Sep 17 00:00:00 2001 From: YL <20972148+LiYunyang@users.noreply.github.com> Date: Mon, 2 Dec 2024 16:10:33 -0600 Subject: [PATCH 4/4] [doc] line width --- maps/python/map_modules.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/maps/python/map_modules.py b/maps/python/map_modules.py index 872ae7d3..fbd72715 100644 --- a/maps/python/map_modules.py +++ b/maps/python/map_modules.py @@ -876,8 +876,8 @@ class ReprojectMaps(object): defined by the mask. If the mask is not provided, it will be determined from the non-zero pixels of the first reprojected map. mask : G3SkyMapMask, G3SkyMap, or np.ndarray, Optional. - Mask to be used for partial reproject. This should be of the same size as the output map. - For numpy array, all zeros/inf/nan/hp.UNSEEN pixels are skipped. + Mask to be used for partial reproject. This should be of the same size as the + output map. For numpy array, all zeros/inf/nan/hp.UNSEEN pixels are skipped. """ def __init__(self, map_stub=None, rebin=1, interp=False, weighted=True,