diff --git a/maps/CMakeLists.txt b/maps/CMakeLists.txt index d0d452c9..1edd89d9 100644 --- a/maps/CMakeLists.txt +++ b/maps/CMakeLists.txt @@ -1,6 +1,7 @@ add_spt3g_library(maps SHARED src/chealpix.c src/G3SkyMap.cxx + src/G3SkyMapMask.cxx src/HealpixSkyMapInfo.cxx src/HealpixSkyMap.cxx src/FlatSkyProjection.cxx @@ -27,6 +28,7 @@ add_spt3g_test(flatsky_maps) add_spt3g_test(flatskymap_operators) add_spt3g_test(healpixskymap_operators) add_spt3g_test(weights_operators) +add_spt3g_test(mask_operators) add_spt3g_test(healpix_maps) add_spt3g_test(fitsio) add_spt3g_test(map_modules) diff --git a/maps/README.rst b/maps/README.rst index 0c4ae86c..bff5b23b 100644 --- a/maps/README.rst +++ b/maps/README.rst @@ -56,7 +56,34 @@ Sparsity By default, both Healpix and flat-sky maps are initialized in sparse mode. This imposes a slight performance penalty but will result in the map storing only non-zero portions (with caveats, see details above), substantially reducing RAM usage. Some map operations, in particular casting to numpy arrays, will result in the implicit conversion of the map to dense storage, which can result in sudden increases in RAM usage. The current sparsity mode can be examined or changed with the ``sparse`` property (flat sky maps) or the ``dense``, ``ringsparse``, or ``indexedsparse`` properties (Healpix maps). Serialization to ``.g3`` files will maintain the current sparsity scheme, as do arithmetic operators where possible. Serialization to ``.fits`` files implicitly converts flat sky maps to dense mode, but preserves the sparsity of Healpix maps. The current number of stored pixels can be obtained using the ``npix_allocated`` property, and the number of non-zero pixels can be obtained using the ``npix_nonzero`` property. Dense maps can be efficiently compactified in memory using the ``G3SkyMap.compact`` method, or the ``CompactMaps`` pipeline module. -Beyond paying attention to implicit conversions to dense storage and the performance impact of sparse storage (which is small), users of this code do not need to worry about the storage mode--all interfaces are identical in all modes. +Beyond paying attention to implicit conversions to dense storage and the performance impact of sparse storage (which is small, at least for FlatSkyMap objects), users of this code do not need to worry about the storage mode--all interfaces are identical in all modes. + +Masking +======= + +Maps containing only boolean data for each pixel are stored as ``G3SkyMapMask`` objects. Such mask objects have a ``.parent`` attribute which is a shallow clone of the map object with which they are associated (to check for shape compatibility). + +Masks are returned when using comparison operators with map objects, e.g. ``map1 > 5`` or ``map1 == map2``. The supported comparison operators are: ``>, >=, ==, !=, <=, <``. Masks can also be combined together using logical operators, e.g. ``mask3 = mask1 & mask2`` or ``mask1 ^= mask2``. The supported comparison operators are: ``&, &=, |, |=, ^, ^=``. Masks can also be checked for equality to other masks using ``==`` and ``!=`` operators. + +Mask objects can be ``clone``'ed the same way as maps. A map can be converted to a boolean mask using ``G3SkyMap.to_mask()``, which returns a mask which is ``True`` wherever the map is non-zero (optionally excluding nan or inf pixels). A mask can be converted back to a map object using ``G3SkyMapMask.to_map()``, which returns a sparse, unit-less, unweighted, unpolarized map object of the same type as ``G3SkyMapMask.parent``, containing double ``1.0`` wherever the mask is ``True``. + +Masks can also be applied to maps or masks using the appropriate ``.apply_mask`` method, with optional inversion; alternatively maps can also be directly multiplied by a compatible mask object. A list of non-zero pixels can be returned using ``.nonzero()`` (note that this returns a single vector of pixel positions), and mask contents can be checked using ``.all()``, ``.any()`` and ``.sum()``. Mask contents can be inverted in-place using ``.invert()``. + +Mask objects cannot be accessed using ``numpy`` slicing, or converted directly to arrays, because ``numpy`` does not represent boolean values as single bits. To be able to use ``numpy`` tools with masks, you need to first convert the mask to a dense map using ``.to_map()``. All associated methods of the parent map are accessible as attributes of the mask object in python, e.g. ``mask.angles_to_pixels()`` works as one would expect. + +Mask Memory Usage +----------------- + +The current implementation of masks is to use a dense ``std::vector`` as the data storage backend, which uses 64x less memory than a dense map (``std::vector``) of the same dimensions. This implementation is sufficient for ``FlatSkyMap`` objects, since these are typically O(50\%) full populated in their sparse state; however, the memory savings for ``HealpixSkyMap`` objects is not as significant when observing sufficiently small patches of sky. Future work would enable a similar sparse storage backend for masks. + +In general, when working with high-resolution maps of any sort, it is important to think carefully about doing the sorts of operations that can balloon memory usage, e.g. taking care to preserve the sparsity of maps by avoiding numpy operations if possible, or using in-place operations to avoid unintentionally creating extra maps or masks in memory. + +Statistics +========== + +Most ``numpy.ufunc``-like methods are defined for map objects, namely ``.all(), .any(), .sum(), .mean(), .median(), .var(), .std(), .min(), .max(), .argmin(), .argmax()``. All methods take an optional ``where`` argument, which can be a compatible ``G3SkyMapMask`` object, or size-compatible 1-D ``numpy`` array that can be converted into one. In addition, these methods are called under the hood when using the numpy equivalent functions (``numpy.all()``, etc), in order to preserve the sparsity of the input map. Methods that ignore ``NaN`` values are also defined (``.nansum()``, etc), which behave much like the standard methods, except that calling ``numpy.nansum()`` and friends on a map object does *not* preserve sparsity. + +Map values can be tested using ``.isnan(), .isinf(), .isfinite()`` methods as well; these return ``G3SkyMapMask`` objects. Map Interpolation ================= diff --git a/maps/include/maps/FlatSkyMap.h b/maps/include/maps/FlatSkyMap.h index fc18ef4e..b75c8503 100644 --- a/maps/include/maps/FlatSkyMap.h +++ b/maps/include/maps/FlatSkyMap.h @@ -11,7 +11,7 @@ #include class DenseMapData; -class SparseMapData; +template class SparseMapData; class FlatSkyMap : public G3FrameObject, public G3SkyMap { public: @@ -73,6 +73,7 @@ class FlatSkyMap : public G3FrameObject, public G3SkyMap { // * virtual G3SkyMap &operator*=(const G3SkyMap &rhs) override; + virtual G3SkyMap &operator*=(const G3SkyMapMask &rhs) override; virtual G3SkyMap &operator*=(double rhs) override; // / @@ -91,6 +92,7 @@ class FlatSkyMap : public G3FrameObject, public G3SkyMap { bool IsCompatible(const G3SkyMap & other) const override; void NonZeroPixels(std::vector &indices, std::vector &data) const; + void ApplyMask(const G3SkyMapMask &mask, bool inverse=false) override; void SetProj(MapProjection proj); void SetAlphaCenter(double alpha); @@ -185,7 +187,7 @@ class FlatSkyMap : public G3FrameObject, public G3SkyMap { FlatSkyProjection proj_info; // projection parameters and functions DenseMapData *dense_; - SparseMapData *sparse_; + SparseMapData *sparse_; uint64_t xpix_, ypix_; bool flat_pol_; diff --git a/maps/include/maps/G3SkyMap.h b/maps/include/maps/G3SkyMap.h index e63ab240..d09a6232 100644 --- a/maps/include/maps/G3SkyMap.h +++ b/maps/include/maps/G3SkyMap.h @@ -7,6 +7,8 @@ #include #include +#include + #include enum MapCoordReference { @@ -97,12 +99,61 @@ class G3SkyMap { // * virtual G3SkyMap &operator*=(const G3SkyMap &rhs); + virtual G3SkyMap &operator*=(const G3SkyMapMask &rhs); virtual G3SkyMap &operator*=(double rhs); // / virtual G3SkyMap &operator/=(const G3SkyMap &rhs); virtual G3SkyMap &operator/=(double rhs); + // Comparison operations + + // < + virtual G3SkyMapMask operator<(const G3SkyMap &rhs); + virtual G3SkyMapMask operator<=(const G3SkyMap &rhs); + virtual G3SkyMapMask operator<(double rhs); + virtual G3SkyMapMask operator<=(double rhs); + + // == + virtual G3SkyMapMask operator==(const G3SkyMap &rhs); + virtual G3SkyMapMask operator!=(const G3SkyMap &rhs); + virtual G3SkyMapMask operator==(double rhs); + virtual G3SkyMapMask operator!=(double rhs); + + // < + virtual G3SkyMapMask operator>(const G3SkyMap &rhs); + virtual G3SkyMapMask operator>=(const G3SkyMap &rhs); + virtual G3SkyMapMask operator>(double rhs); + virtual G3SkyMapMask operator>=(double rhs); + + virtual bool all(G3SkyMapMaskConstPtr where=NULL) const; + virtual bool any(G3SkyMapMaskConstPtr where=NULL) const; + virtual double sum(G3SkyMapMaskConstPtr where=NULL) const; + virtual double mean(G3SkyMapMaskConstPtr where=NULL) const; + virtual double median(G3SkyMapMaskConstPtr where=NULL) const; + virtual double var(size_t ddof, G3SkyMapMaskConstPtr where=NULL) const; + virtual double std(size_t ddof, G3SkyMapMaskConstPtr where=NULL) const; + virtual double min(G3SkyMapMaskConstPtr where=NULL) const; + virtual double max(G3SkyMapMaskConstPtr where=NULL) const; + virtual size_t argmin(G3SkyMapMaskConstPtr where=NULL) const; + virtual size_t argmax(G3SkyMapMaskConstPtr where=NULL) const; + + virtual double nansum(G3SkyMapMaskConstPtr where=NULL) const; + virtual double nanmean(G3SkyMapMaskConstPtr where=NULL) const; + virtual double nanmedian(G3SkyMapMaskConstPtr where=NULL) const; + virtual double nanvar(size_t ddof, G3SkyMapMaskConstPtr where=NULL) const; + virtual double nanstd(size_t ddof, G3SkyMapMaskConstPtr where=NULL) const; + virtual double nanmin(G3SkyMapMaskConstPtr where=NULL) const; + virtual double nanmax(G3SkyMapMaskConstPtr where=NULL) const; + virtual size_t nanargmin(G3SkyMapMaskConstPtr where=NULL) const; + virtual size_t nanargmax(G3SkyMapMaskConstPtr where=NULL) const; + + virtual G3SkyMapMask isinf(G3SkyMapMaskConstPtr where=NULL) const; + virtual G3SkyMapMask isnan(G3SkyMapMaskConstPtr where=NULL) const; + virtual G3SkyMapMask isfinite(G3SkyMapMaskConstPtr where=NULL) const; + + virtual void ApplyMask(const G3SkyMapMask &mask, bool inverse=false); + // Pointing information std::vector AnglesToPixels(const std::vector & alphas, const std::vector & deltas) const; @@ -150,6 +201,9 @@ class G3SkyMap { throw std::runtime_error("Compactification not implemented"); } + virtual G3SkyMapMaskPtr MakeMask(bool zero_nans = false, + bool zero_infs = false) const; + protected: MapPolConv pol_conv_; @@ -362,6 +416,7 @@ class G3SkyMapWeights : public G3FrameObject { // Mask G3SkyMapWeights &operator*=(const G3SkyMap &rhs); + G3SkyMapWeights &operator*=(const G3SkyMapMask &rhs); G3SkyMapPtr Det() const; G3SkyMapPtr Cond() const; @@ -369,6 +424,8 @@ class G3SkyMapWeights : public G3FrameObject { boost::shared_ptr Rebin(size_t scale) const; + void ApplyMask(const G3SkyMapMask &mask, bool inverse=false); + void Compact(bool zero_nans = false); boost::shared_ptr Clone(bool copy_data) const { diff --git a/maps/include/maps/G3SkyMapMask.h b/maps/include/maps/G3SkyMapMask.h new file mode 100644 index 00000000..7d56fe43 --- /dev/null +++ b/maps/include/maps/G3SkyMapMask.h @@ -0,0 +1,128 @@ +#ifndef _G3_SKYMAPMASK_H +#define _G3_SKYMAPMASK_H + +#include + +#include + +/* + * The following implements a companion object to a G3SkyMap containing + * booleans for whether to use (true) or ignore (false) certain pixels. + * The pixelization is stored externally in the parent map, as is the + * (potential) translation from 2D to 1D pixel indices. + * + * By default, initializes to all zeroes. If use_data = True, it will + * interpret the input sky map both as projection information and as a + * source of data, initializing the mask to True where the input map is + * non-zero. + */ + +class G3SkyMap; + +class G3SkyMapMask : public G3FrameObject { +public: + G3SkyMapMask(const G3SkyMap &parent, bool use_data = false, + bool zero_nans = false, bool zero_infs = false); + G3SkyMapMask(const G3SkyMap &parent, boost::python::object v, + bool zero_nans = false, bool zero_infs = false); + G3SkyMapMask(const G3SkyMapMask &); + virtual ~G3SkyMapMask() {}; + + // Copy + boost::shared_ptr Clone(bool copy_data = true) const; + boost::shared_ptr ArrayClone(boost::python::object v, + bool zero_nans = false, bool zero_infs = false) const; + + // Return a (modifiable) pixel value + std::vector::reference operator [] (size_t i); + bool operator [] (size_t i) const { return this->at(i); }; + bool at(size_t i) const; + + // Logic operations: + + G3SkyMapMask &operator |=(const G3SkyMapMask &rhs); + G3SkyMapMask &operator &=(const G3SkyMapMask &rhs); + G3SkyMapMask &operator ^=(const G3SkyMapMask &rhs); + G3SkyMapMask &invert(); // Basically ~= + bool all() const; + bool any() const; + size_t sum() const; + size_t size() const; + std::vector NonZeroPixels() const; + + G3SkyMapMask operator ~() const; + G3SkyMapMask operator |(const G3SkyMapMask &rhs) const; + G3SkyMapMask operator &(const G3SkyMapMask &rhs) const; + G3SkyMapMask operator ^(const G3SkyMapMask &rhs) const; + G3SkyMapMask operator ==(const G3SkyMapMask &rhs) const; + G3SkyMapMask operator !=(const G3SkyMapMask &rhs) const; + + void ApplyMask(const G3SkyMapMask &rhs, bool inverse=false); + + // Information + bool IsCompatible(const G3SkyMap &map) const; + bool IsCompatible(const G3SkyMapMask &mask) const; + bool IsDense() const { return true; } + + // The map for projection info + boost::shared_ptr Parent() const { return parent_; } + + // Return a 1-or-0 sky-map with the contents of the mask (e.g. for plotting) + boost::shared_ptr MakeBinaryMap() const; + + class const_iterator { + public: + typedef std::pair value_type; + typedef value_type & reference; + typedef value_type * pointer; + + const_iterator(const G3SkyMapMask &mask, bool begin); + + bool operator==(const const_iterator & other) const { + return index_ == other.index_; + }; + + bool operator!=(const const_iterator & other) const { + return index_ != other.index_; + }; + + reference operator*() { return value_; }; + pointer operator->() { return &value_; }; + + const_iterator operator++(); + const_iterator operator++(int) { const_iterator i = *this; ++(*this); return i; } + + private: + uint64_t index_; + value_type value_; + const G3SkyMapMask &mask_; + + void set_value() { + value_.first = index_; + value_.second = mask_.at(index_); + } + }; + + const_iterator begin() const { return const_iterator(*this, true); }; + const_iterator end() const { return const_iterator(*this, false); }; + +private: + G3SkyMapMask() {} // Fake out for serialization + template void serialize(A &ar, const unsigned v); + friend class cereal::access; + + std::vector data_; + boost::shared_ptr parent_; + + // Populate + void FillFromMap(const G3SkyMap &m, bool zero_nans = false, bool zero_infs = false); + void FillFromArray(boost::python::object v, bool zero_nans = false, bool zero_infs = false); + + SET_LOGGER("G3SkyMapMask"); +}; + +G3_POINTERS(G3SkyMapMask); +G3_SERIALIZABLE(G3SkyMapMask, 1); + +#endif + diff --git a/maps/include/maps/HealpixSkyMap.h b/maps/include/maps/HealpixSkyMap.h index a96f56d8..e306dea6 100644 --- a/maps/include/maps/HealpixSkyMap.h +++ b/maps/include/maps/HealpixSkyMap.h @@ -10,7 +10,7 @@ #include #include -class SparseMapData; +template class SparseMapData; class HealpixSkyMap : public G3FrameObject, public G3SkyMap { @@ -53,6 +53,7 @@ class HealpixSkyMap : public G3FrameObject, public G3SkyMap { // * virtual G3SkyMap &operator*=(const G3SkyMap &rhs) override; + virtual G3SkyMap &operator*=(const G3SkyMapMask &rhs) override; virtual G3SkyMap &operator*=(double rhs) override; // / @@ -71,6 +72,7 @@ class HealpixSkyMap : public G3FrameObject, public G3SkyMap { bool IsCompatible(const G3SkyMap & other) const override; void NonZeroPixels(std::vector &indices, std::vector &data) const; // Iterators better? + void ApplyMask(const G3SkyMapMask &mask, bool inverse=false) override; size_t nside() const {return info_.nside();} bool nested() const {return info_.nested();} @@ -141,7 +143,7 @@ class HealpixSkyMap : public G3FrameObject, public G3SkyMap { private: HealpixSkyMapInfo info_; std::vector *dense_; - SparseMapData *ring_sparse_; + SparseMapData *ring_sparse_; std::unordered_map *indexed_sparse_; SET_LOGGER("HealpixSkyMap"); diff --git a/maps/include/maps/maputils.h b/maps/include/maps/maputils.h index ee93ce52..c8935ca8 100644 --- a/maps/include/maps/maputils.h +++ b/maps/include/maps/maputils.h @@ -8,14 +8,6 @@ #include #include -// Return a map that is set to one where the input map is non-zero. -// Optionally zero out nans and/or infs if zero_nans or zero_infs is true. -G3SkyMapPtr GetMaskMap(G3SkyMapConstPtr m, bool zero_nans=false, bool zero_infs=false); - -// Zero out any values in the input map where the mask is zero -// If inverse is true, zero the input map where the mask is non-zero -void ApplyMask(G3SkyMapPtr m, G3SkyMapConstPtr mask, bool inverse=false); - // Divide out map weights from the sky maps void RemoveWeightsT(G3SkyMapPtr T, G3SkyMapWeightsConstPtr W, bool zero_nans = false); void RemoveWeights(G3SkyMapPtr T, G3SkyMapPtr Q, G3SkyMapPtr U, G3SkyMapWeightsConstPtr W, @@ -30,7 +22,7 @@ boost::python::tuple GetRaDecMap(G3SkyMapConstPtr m); // Return a mask that is one for all pixels in the input map that are within the // given ra and dec rectangular bounds, and zero otherwise -G3SkyMapPtr GetRaDecMask(G3SkyMapConstPtr m, double ra_left, double ra_right, +G3SkyMapMaskPtr GetRaDecMask(G3SkyMapConstPtr m, double ra_left, double ra_right, double dec_bottom, double dec_top); // Reproject the input map onto the output map grid, optionally @@ -44,30 +36,22 @@ void FlattenPol(FlatSkyMapPtr Q, FlatSkyMapPtr U, G3SkyMapWeightsPtr W=NULL, dou // Compute map moments up to fourth order (mean, var, skewness, kurtosis) of the input map, // optionally excluding any pixels that are zero in the input mask, or zero/nan/inf in the input map -std::vector GetMapStats(G3SkyMapConstPtr m, G3SkyMapConstPtr mask=NULL, int order=2, - bool ignore_zeros=false, bool ignore_nans=false, bool ignore_infs=false); - -// Compute the median of the input map, optionally excluding any pixels that are -// zero in the input mask, or zero/nan/inf in the input map -double GetMapMedian(G3SkyMapConstPtr m, G3SkyMapConstPtr mask=NULL, - bool ignore_zeros=false, bool ignore_nans=false, bool ignore_infs=false); - -// Compute the min and max of the values in the input map, optionally excluding -// any pixels that are zero in the input mask, or zero/nan/inf in the input map -std::vector GetMapMinMax(G3SkyMapConstPtr m, G3SkyMapConstPtr mask=NULL, +std::vector GetMapMoments(G3SkyMapConstPtr m, G3SkyMapMaskConstPtr mask=NULL, int order=2, bool ignore_zeros=false, bool ignore_nans=false, bool ignore_infs=false); // Compute the histogram of the input map pixels, grouping the values into bins // defined by the array of bin edges, and optionally excluding any pixels that are // zero in the input mask, or zero/nan/inf in the input map std::vector GetMapHist(G3SkyMapConstPtr m, const std::vector &bin_edges, - G3SkyMapConstPtr mask, bool ignore_zeros, bool ignore_nans, bool ignore_infs); + G3SkyMapMaskConstPtr mask=NULL, bool ignore_zeros=false, bool ignore_nans=false, + bool ignore_infs=false); // Convolve the input flat sky map with a filter kernel FlatSkyMapPtr ConvolveMap(FlatSkyMapConstPtr map, FlatSkyMapConstPtr kernel); // Point source masking -void MakePointSourceMask(G3SkyMapPtr map, const std::vector & ra, +G3SkyMapMaskPtr +MakePointSourceMask(G3SkyMapConstPtr map, const std::vector & ra, const std::vector & dec, const std::vector & radius); // Python bindings diff --git a/maps/python/map_modules.py b/maps/python/map_modules.py index a4f429dc..d4d2c8a6 100644 --- a/maps/python/map_modules.py +++ b/maps/python/map_modules.py @@ -158,7 +158,7 @@ def MakeMapsPolarized(frame, pol_conv=maps.MapPolConv.IAU): umap.pol_type = maps.MapPolType.U umap.pol_conv = pol_conv frame["U"] = umap - mask = maps.get_mask_map(wgt) + mask = wgt.to_mask().to_map() wgt_out = maps.G3SkyMapWeights(frame["T"], polarized=True) wgt_out.TT = wgt diff --git a/maps/python/skymapaddons.py b/maps/python/skymapaddons.py index 2b595fdf..2f157232 100644 --- a/maps/python/skymapaddons.py +++ b/maps/python/skymapaddons.py @@ -1,6 +1,6 @@ import numpy import warnings -from spt3g.maps import G3SkyMapWeights, G3SkyMap, FlatSkyMap +from spt3g.maps import G3SkyMapWeights, G3SkyMap, FlatSkyMap, G3SkyMapMask # This file adds extra functionality to the python interface to G3SkyMap and # G3SkyMapWeights. This is done in ways that exploit a large fraction of @@ -29,6 +29,31 @@ def numpycompat(a, b, op): for x in ['__and__', '__divmod__', '__floordiv__', '__iand__', '__ifloordiv__', '__imod__', '__ior__', '__mod__', '__or__', '__rdivmod__', '__rmod__', '__rpow__']: setattr(G3SkyMap, x, lambda a, b, op=x: numpycompat(a, b, op)) +for op in ["all", "any", "sum", "mean", "var", "std", "min", "max", "argmin", "argmax"]: + def ufunc_wrapper(op): + def ufunc(a, *args, **kwargs): + bound_args = {} + bound = getattr(a, "_c" + op) + if op in ["std", "var", "nanstd", "nanvar"]: + if "ddof" in kwargs: + bound_args["ddof"] = kwargs.pop("ddof") + if "where" in kwargs: + where = kwargs.pop("where") + if isinstance(where, numpy.ndarray): + where = G3SkyMapMask(a, where) + bound_args["where"] = where + for k in ["axis", "out", "dtype"]: + if k in kwargs and kwargs[k] is None: + kwargs.pop(k) + if len(args) or len(kwargs): + raise TypeError("ufunc {} does not support complex arguments".format(op)) + return bound(**bound_args) + ufunc.__doc__ = getattr(numpy.ndarray, op).__doc__ + return ufunc + setattr(G3SkyMap, op, ufunc_wrapper(op)) + if op in ["all", "any", "sum"]: + setattr(G3SkyMapMask, op, ufunc_wrapper(op)) + # Make weight maps so that you can index them and get the full 3x3 weight matrix def skymapweights_keys(self): @@ -152,7 +177,9 @@ def skymapweights_reshape(self, width, height, fill=0): def skymapweights_getattr(self, x): if x == "flat_pol" and isinstance(self.TQ, FlatSkyMap): return getattr(self.TQ, x) - return getattr(self.TT, x) + if hasattr(self.TT, x): + return getattr(self.TT, x) + raise AttributeError("'{}' object has no attribute '{}'".format(type(self), x)) G3SkyMapWeights.__getattr__ = skymapweights_getattr del skymapweights_getattr oldsetattr = G3SkyMapWeights.__setattr__ @@ -160,16 +187,27 @@ def skymapweights_setattr(self, x, val): try: oldsetattr(self, x, val) except AttributeError: - setattr(self.TT, x, val) - if self.polarized: + if self.TT is not None: + setattr(self.TT, x, val) + if self.TQ is not None: setattr(self.TQ, x, val) + if self.TU is not None: setattr(self.TU, x, val) + if self.QQ is not None: setattr(self.QQ, x, val) + if self.QU is not None: setattr(self.QU, x, val) + if self.UU is not None: setattr(self.UU, x, val) G3SkyMapWeights.__setattr__ = skymapweights_setattr del skymapweights_setattr +def skymapmask_getattr(self, x): + if hasattr(self.parent, x): + return getattr(self.parent, x) + raise AttributeError("'{}' object has no attribute '{}'".format(type(self), x)) +G3SkyMapMask.__getattr__ = skymapmask_getattr +del skymapmask_getattr def skymap_clone(self, copy_data=True): with warnings.catch_warnings(): diff --git a/maps/src/FlatSkyMap.cxx b/maps/src/FlatSkyMap.cxx index 478253d7..f56c98e1 100644 --- a/maps/src/FlatSkyMap.cxx +++ b/maps/src/FlatSkyMap.cxx @@ -8,6 +8,7 @@ #include #include +#include #include "mapdata.h" @@ -81,7 +82,7 @@ FlatSkyMap::FlatSkyMap(const FlatSkyMap & fm) : if (fm.dense_) dense_ = new DenseMapData(*fm.dense_); else if (fm.sparse_) - sparse_ = new SparseMapData(*fm.sparse_); + sparse_ = new SparseMapData(*fm.sparse_); } FlatSkyMap::~FlatSkyMap() @@ -159,7 +160,7 @@ FlatSkyMap::load(A &ar, unsigned v) assert(dense_->ydim() == ypix_); break; case 1: - sparse_ = new SparseMapData(0, 0); + sparse_ = new SparseMapData(0, 0); ar & make_nvp("sparse", *sparse_); assert(sparse_->xdim() == xpix_); assert(sparse_->ydim() == ypix_); @@ -285,7 +286,7 @@ FlatSkyMap::const_iterator::operator++() x_ = iter.x; y_ = iter.y; } else if (map_.sparse_) { - SparseMapData::const_iterator iter(*map_.sparse_, x_, y_); + SparseMapData::const_iterator iter(*map_.sparse_, x_, y_); ++iter; x_ = iter.x; y_ = iter.y; @@ -318,7 +319,7 @@ FlatSkyMap::ConvertToSparse() if (!dense_) return; - sparse_ = new SparseMapData(*dense_); + sparse_ = new SparseMapData(*dense_); delete dense_; dense_ = NULL; } @@ -369,7 +370,7 @@ FlatSkyMap::operator () (size_t x, size_t y) if (dense_) return (*dense_)(x, y); if (!sparse_) - sparse_ = new SparseMapData(xpix_, ypix_); + sparse_ = new SparseMapData(xpix_, ypix_); return (*sparse_)(x, y); } @@ -413,7 +414,7 @@ G3SkyMap &FlatSkyMap::operator op(const G3SkyMap &rhs) {\ ConvertToDense(); \ (*dense_) op (*b.dense_); \ } else if (b.sparse_) { \ - sparse_ = new SparseMapData(xpix_, ypix_); \ + sparse_ = new SparseMapData(xpix_, ypix_); \ (*sparse_) op (*b.sparse_); \ } else { \ sparsenull \ @@ -473,6 +474,19 @@ G3SkyMap &FlatSkyMap::operator *=(const G3SkyMap &rhs) { } } +G3SkyMap & +FlatSkyMap::operator *=(const G3SkyMapMask &rhs) +{ + g3_assert(rhs.IsCompatible(*this)); + + for (auto i : *this) { + if (!rhs.at(i.first) && i.second != 0) + (*this)[i.first] = 0; + } + + return *this; +} + G3SkyMap & FlatSkyMap::operator +=(double b) { @@ -625,6 +639,16 @@ FlatSkyMap::NonZeroPixels(std::vector &indices, } } +void +FlatSkyMap::ApplyMask(const G3SkyMapMask &mask, bool inverse) +{ + g3_assert(mask.IsCompatible(*this)); + + for (auto i: *this) + if (i.second != 0 && mask.at(i.first) == inverse) + (*this)[i.first] = 0; +} + std::vector FlatSkyMap::shape() const { return {xpix_, ypix_}; } @@ -1017,6 +1041,45 @@ flatskymap_setitem_1d(G3SkyMap &skymap, size_t i, double val) skymap[i] = val; } +static std::vector +flatskymap_getitem_masked(const FlatSkyMap &skymap, const G3SkyMapMask &m) +{ + g3_assert(m.IsCompatible(skymap)); + std::vector out; + + for (auto i : skymap) { + if (m.at(i.first)) + out.push_back(i.second); + } + + return out; +} + +static void +flatskymap_setitem_masked(FlatSkyMap &skymap, const G3SkyMapMask &m, + bp::object val) +{ + g3_assert(m.IsCompatible(skymap)); + + if (bp::extract(val).check()) { + double dval = bp::extract(val)(); + for (auto i : skymap) { + if (m.at(i.first)) + skymap[i.first] = dval; + } + } else { + // XXX: the iterable case probably be optimized for numpy arrays + // XXX: check for size congruence first? + size_t j = 0; + for (auto i : skymap) { + if (m.at(i.first)) { + skymap[i.first] = bp::extract(val[j])(); + j++; + } + } + } +} + static bool flatskymap_pysparsity_get(const FlatSkyMap &fsm) { @@ -1214,6 +1277,8 @@ PYBINDINGS("maps") .def("__getitem__", flatskymap_getitem_2d) .def("__setitem__", flatskymap_setitem_2d) .def("__setitem__", flatskymap_setslice_1d) + .def("__getitem__", flatskymap_getitem_masked) + .def("__setitem__", flatskymap_setitem_masked) ; register_pointer_conversions(); @@ -1226,8 +1291,6 @@ PYBINDINGS("maps") #endif implicitly_convertible(); - implicitly_convertible(); implicitly_convertible(); - implicitly_convertible(); } diff --git a/maps/src/G3SkyMap.cxx b/maps/src/G3SkyMap.cxx index bf7f4cb3..24b668ef 100644 --- a/maps/src/G3SkyMap.cxx +++ b/maps/src/G3SkyMap.cxx @@ -2,9 +2,11 @@ #include #include #include +#include #include #include #include +#include namespace bp=boost::python; @@ -236,6 +238,13 @@ size_t G3SkyMap::size() const return s; } +G3SkyMapMaskPtr +G3SkyMap::MakeMask(bool zero_nans, bool zero_infs) const +{ + G3SkyMapMaskPtr m(new G3SkyMapMask(*this, true, zero_nans, zero_infs)); + return m; +} + G3SkyMap &G3SkyMap::operator+=(const G3SkyMap & rhs) { g3_assert(IsCompatible(rhs)); @@ -285,6 +294,17 @@ G3SkyMap &G3SkyMap::operator*=(const G3SkyMap &rhs) return *this; } +G3SkyMap &G3SkyMap::operator*=(const G3SkyMapMask &rhs) +{ + g3_assert(rhs.IsCompatible(*this)); + + for (size_t i = 0; i < size(); i++) { + if (!rhs.at(i) && this->at(i) != 0) + (*this)[i] = 0; + } + return *this; +} + G3SkyMap &G3SkyMap::operator*=(double rhs) { for (size_t i = 0; i < size(); i++) @@ -520,6 +540,26 @@ skymap_pynoninplace(multd, *=, double) skymap_pynoninplace(div, /=, G3SkyMap &) skymap_pynoninplace(divd, /=, double) +static G3SkyMapPtr +pyskymap_multm(const G3SkyMap &a, const G3SkyMapMask &b) +{ + g3_assert(b.IsCompatible(a)); + G3SkyMapPtr rv = a.Clone(false); + for (auto i: b) { + if (b.at(i.first) && a.at(i.first) != 0) + (*rv)[i.first] = a.at(i.first); + } + + return rv; +} + +static G3SkyMapPtr +pyskymap_imultm(G3SkyMapPtr a, const G3SkyMapMask &b) +{ + (*a) *= b; + return a; +} + static G3SkyMapPtr pyskymap_rsubd(const G3SkyMap &a, const double b) { @@ -597,68 +637,386 @@ pyskymap_pow(const G3SkyMap &a, const G3SkyMap &b) return rv; } -#define skymap_comp(name, oper) \ -static G3SkyMapPtr \ -pyskymap_##name(const G3SkyMap &a, const G3SkyMap &b) \ +#define skymap_comp(oper) \ +G3SkyMapMask \ +G3SkyMap::operator oper(const G3SkyMap &rhs) \ { \ - g3_assert(a.IsCompatible(b)); \ - g3_assert(a.units == b.units); \ - G3SkyMapPtr rv = a.Clone(false); \ - rv->units = G3Timestream::None; \ - rv->weighted = false; \ - for (size_t i = 0; i < a.size(); i++) { \ - if (a.at(i) oper b.at(i)) \ - (*rv)[i] = 1; \ + g3_assert(IsCompatible(rhs)); \ + g3_assert(units == rhs.units); \ + G3SkyMapMask rv(*this); \ + for (size_t i = 0; i < size(); i++) { \ + if (at(i) oper rhs.at(i)) \ + rv[i] = true; \ } \ return rv; \ } -skymap_comp(lt, <) -skymap_comp(le, <=) -skymap_comp(eq, ==) -skymap_comp(ne, !=) -skymap_comp(ge, >=) -skymap_comp(gt, >) +skymap_comp(<) +skymap_comp(<=) +skymap_comp(==) +skymap_comp(!=) +skymap_comp(>=) +skymap_comp(>) -#define skymap_compd(name, oper) \ -static G3SkyMapPtr \ -pyskymap_##name(const G3SkyMap &a, const double b) \ +#define skymap_compd(oper) \ +G3SkyMapMask \ +G3SkyMap::operator oper(double rhs) \ { \ - G3SkyMapPtr rv = a.Clone(false); \ - rv->units = G3Timestream::None; \ - rv->weighted = false; \ - for (size_t i = 0; i < a.size(); i++) { \ - if (a.at(i) oper b) \ - (*rv)[i] = 1; \ + G3SkyMapMask rv(*this); \ + for (size_t i = 0; i < size(); i++) { \ + if (at(i) oper rhs) \ + rv[i] = true; \ } \ return rv; \ } -skymap_compd(ltd, <) -skymap_compd(led, <=) -skymap_compd(eqd, ==) -skymap_compd(ned, !=) -skymap_compd(ged, >=) -skymap_compd(gtd, >) +skymap_compd(<) +skymap_compd(<=) +skymap_compd(==) +skymap_compd(!=) +skymap_compd(>=) +skymap_compd(>) -static bool -pyskymap_all(const G3SkyMap &a) +bool +G3SkyMap::all(G3SkyMapMaskConstPtr where) const { - for (size_t i = 0; i < a.size(); i++) - if (a.at(i) == 0) - return false; + if (!where) { + for (size_t i = 0; i < size(); i++) { + if (at(i) == 0) + return false; + } + } else { + g3_assert(where->IsCompatible(*this)); + + for (size_t i = 0; i < size(); i++) { + if (!where->at(i)) + continue; + if (at(i) == 0) + return false; + } + } + return true; } -static bool -pyskymap_any(const G3SkyMap &a) +bool +G3SkyMap::any(G3SkyMapMaskConstPtr where) const { - for (size_t i = 0; i < a.size(); i++) - if (a.at(i) != 0) - return true; + if (!where) { + for (size_t i = 0; i < size(); i++) { + if (at(i) != 0) + return true; + } + } else { + g3_assert(where->IsCompatible(*this)); + + for (size_t i = 0; i < size(); i++) { + if (!where->at(i)) + continue; + if (at(i) != 0) + return true; + } + } + return false; } +double +G3SkyMap::sum(G3SkyMapMaskConstPtr where) const +{ + double s = 0; + + if (!where) { + for (size_t i = 0; i < size(); i++) { + s += at(i); + } + } else { + g3_assert(where->IsCompatible(*this)); + + for (size_t i = 0; i < size(); i++) { + if (!where->at(i)) + continue; + s += at(i); + } + } + + return s; +} + +double +G3SkyMap::mean(G3SkyMapMaskConstPtr where) const +{ + double m = 0; + size_t n = 0; + + if (!where) { + n = size(); + for (size_t i = 0; i < n; i++) { + m += at(i); + } + } else { + g3_assert(where->IsCompatible(*this)); + + for (size_t i = 0; i < size(); i++) { + if (!where->at(i)) + continue; + m += at(i); + n++; + } + } + + return m / n; +} + +double +G3SkyMap::var(size_t ddof, G3SkyMapMaskConstPtr where) const +{ + double m1 = 0; + double m2 = 0; + double a, b; + size_t n = 0; + + if (!where) { + n = size(); + for (size_t i = 0; i < n; i++) { + a = at(i); + m1 += a; + m2 += a * a; + } + } else { + g3_assert(where->IsCompatible(*this)); + + for (size_t i = 0; i < size(); i++) { + if (!where->at(i)) + continue; + n++; + a = at(i); + m1 += a; + m2 += a * a; + } + } + + return (m2 - m1 * m1 / (double)n) / (double)(n - ddof); +} + +double +G3SkyMap::std(size_t ddof, G3SkyMapMaskConstPtr where) const +{ + return sqrt(var(ddof, where)); +} + +double +G3SkyMap::min(G3SkyMapMaskConstPtr where) const +{ + double m = std::numeric_limits::infinity(); + double v; + + if (!where) { + for (size_t i = 0; i < size(); i++) { + v = at(i); + if (v < m) + m = v; + } + } else { + g3_assert(where->IsCompatible(*this)); + + for (size_t i = 0; i < size(); i++) { + if (!where->at(i)) + continue; + v = at(i); + if (v < m) + m = v; + } + } + + return m; +} + +size_t +G3SkyMap::argmin(G3SkyMapMaskConstPtr where) const +{ + double m = std::numeric_limits::infinity(); + double v; + size_t j = 0; + + if (!where) { + for (size_t i = 0; i < size(); i++) { + v = at(i); + if (v < m) { + m = v; + j = i; + } + } + } else { + g3_assert(where->IsCompatible(*this)); + + for (size_t i = 0; i < size(); i++) { + if (!where->at(i)) + continue; + v = at(i); + if (v < m) { + m = v; + j = i; + } + } + } + + return j; +} + +double +G3SkyMap::max(G3SkyMapMaskConstPtr where) const +{ + double m = -1 * std::numeric_limits::infinity(); + double v; + + if (!where) { + for (size_t i = 0; i < size(); i++) { + v = at(i); + if (v > m) + m = v; + } + } else { + g3_assert(where->IsCompatible(*this)); + + for (size_t i = 0; i < size(); i++) { + if (!where->at(i)) + continue; + v = at(i); + if (v > m) + m = v; + } + } + + return m; +} + +size_t +G3SkyMap::argmax(G3SkyMapMaskConstPtr where) const +{ + double m = -1 * std::numeric_limits::infinity(); + double v; + size_t j = 0; + + if (!where) { + for (size_t i = 0; i < size(); i++) { + v = at(i); + if (v > m) { + m = v; + j = i; + } + } + } else { + g3_assert(where->IsCompatible(*this)); + + for (size_t i = 0; i < size(); i++) { + if (!where->at(i)) + continue; + v = at(i); + if (v > m) { + m = v; + j = i; + } + } + } + + return j; +} + +double +G3SkyMap::median(G3SkyMapMaskConstPtr where) const +{ + std::vector data; + size_t npix = where ? where->sum() : size(); + if (npix == 0) + return 0; + + if (!where) { + for (size_t i = 0; i < size(); i++) { + data.push_back(at(i)); + } + } else { + g3_assert(where->IsCompatible(*this)); + + for (size_t i = 0; i < size(); i++) { + if (!where->at(i)) + continue; + data.push_back(at(i)); + } + } + + npix = data.size(); + size_t n = npix / 2; + + std::nth_element(data.begin(), data.begin() + n, data.end()); + + // odd-length array + if (npix % 2) + return data[n]; + + // even-length array + double m = data[n]; + std::nth_element(data.begin(), data.begin() + n - 1, data.end()); + return (m + data[n - 1]) / 2.; +} + +#define skymap_test(oper) \ +G3SkyMapMask \ +G3SkyMap::oper(G3SkyMapMaskConstPtr where) const \ +{ \ + G3SkyMapMask mask(*this); \ + if (!where) { \ + for (size_t i = 0; i < size(); i++) { \ + if (std::oper(at(i))) \ + mask[i] = true; \ + } \ + } else { \ + g3_assert(where->IsCompatible(*this)); \ + for (size_t i = 0; i < size(); i++) { \ + if (!where->at(i)) \ + continue; \ + if (std::oper(at(i))) \ + mask[i] = true; \ + } \ + } \ + return mask; \ +} + +skymap_test(isinf); +skymap_test(isnan); +skymap_test(isfinite); + +#define skymap_nanoper(oper, type) \ +type \ +G3SkyMap::nan##oper(G3SkyMapMaskConstPtr where) const \ +{ \ + G3SkyMapMask mask = isnan(where); \ + mask.invert(); \ + return oper(boost::make_shared(mask)); \ +} + +skymap_nanoper(sum, double); +skymap_nanoper(mean, double); +skymap_nanoper(median, double); +skymap_nanoper(min, double); +skymap_nanoper(max, double); +skymap_nanoper(argmin, size_t); +skymap_nanoper(argmax, size_t); + +double +G3SkyMap::nanvar(size_t ddof, G3SkyMapMaskConstPtr where) const +{ + G3SkyMapMask mask = isnan(where); + mask.invert(); + return var(ddof, boost::make_shared(mask)); +} + +double +G3SkyMap::nanstd(size_t ddof, G3SkyMapMaskConstPtr where) const +{ + return sqrt(nanvar(ddof, where)); +} + static bool pyskymap_bool(G3SkyMap &skymap) { @@ -669,20 +1027,59 @@ pyskymap_bool(G3SkyMap &skymap) return false; } +static std::vector +pyskymap_nonzero(G3SkyMap &skymap) +{ + return skymap.MakeMask()->NonZeroPixels(); +} + +void +G3SkyMap::ApplyMask(const G3SkyMapMask &mask, bool inverse) +{ + g3_assert(mask.IsCompatible(*this)); + + for (size_t i = 0; i < size(); i++) { + if (!at(i)) + continue; + if (mask.at(i) == inverse) + (*this)[i] = 0; + } +} + +void +G3SkyMapWeights::ApplyMask(const G3SkyMapMask &mask, bool inverse) +{ + if (!!TT) + TT->ApplyMask(mask, inverse); + if (!!TQ) + TQ->ApplyMask(mask, inverse); + if (!!TU) + TU->ApplyMask(mask, inverse); + if (!!QQ) + QQ->ApplyMask(mask, inverse); + if (!!QU) + QU->ApplyMask(mask, inverse); + if (!!UU) + UU->ApplyMask(mask, inverse); +} + G3SkyMapWeights & G3SkyMapWeights::operator+=(const G3SkyMapWeights &rhs) { g3_assert(IsPolarized() == rhs.IsPolarized()); - *TT += *(rhs.TT); - - if (IsPolarized()) { + if (!!TT) + *TT += *rhs.TT; + if (!!TQ) *TQ += *rhs.TQ; + if (!!TU) *TU += *rhs.TU; + if (!!QQ) *QQ += *rhs.QQ; + if (!!QU) *QU += *rhs.QU; + if (!!UU) *UU += *rhs.UU; - } return *this; } @@ -690,18 +1087,23 @@ G3SkyMapWeights::operator+=(const G3SkyMapWeights &rhs) #define skymapweights_inplace(op, rhs_type) \ G3SkyMapWeights &G3SkyMapWeights::operator op(rhs_type rhs) \ { \ - *TT op rhs; \ - if (IsPolarized()) { \ + if (!!TT) \ + *TT op rhs; \ + if (!!TQ) \ *TQ op rhs; \ + if (!!TU) \ *TU op rhs; \ + if (!!QQ) \ *QQ op rhs; \ + if (!!QU) \ *QU op rhs; \ + if (!!UU) \ *UU op rhs; \ - } \ return *this; \ } skymapweights_inplace(*=, const G3SkyMap &); +skymapweights_inplace(*=, const G3SkyMapMask &); skymapweights_inplace(*=, double); skymapweights_inplace(/=, double); @@ -719,6 +1121,45 @@ skymapweights_pynoninplace(multm, *=, G3SkyMap &); skymapweights_pynoninplace(multd, *=, double); skymapweights_pynoninplace(divd, /=, double); +static G3SkyMapWeightsPtr +pyskymapweights_multma(const G3SkyMapWeights &a, const G3SkyMapMask &b) +{ + G3SkyMapWeightsPtr rv(new G3SkyMapWeights()); + if (!!a.TT) + rv->TT = pyskymap_multm(*a.TT, b); + if (!!a.TQ) + rv->TQ = pyskymap_multm(*a.TQ, b); + if (!!a.TU) + rv->TU = pyskymap_multm(*a.TU, b); + if (!!a.QQ) + rv->QQ = pyskymap_multm(*a.QQ, b); + if (!!a.QU) + rv->QU = pyskymap_multm(*a.QU, b); + if (!!a.UU) + rv->UU = pyskymap_multm(*a.UU, b); + + return rv; +} + +static G3SkyMapWeightsPtr +pyskymapweights_imultma(G3SkyMapWeightsPtr a, const G3SkyMapMask &b) +{ + if (!!a->TT) + (*a->TT) *= b; + if (!!a->TQ) + (*a->TQ) *= b; + if (!!a->TU) + (*a->TU) *= b; + if (!!a->QQ) + (*a->QQ) *= b; + if (!!a->QU) + (*a->QU) *= b; + if (!!a->UU) + (*a->UU) *= b; + + return a; +} + MuellerMatrix MuellerMatrix::Inv() const { MuellerMatrix m; @@ -837,13 +1278,16 @@ G3SkyMapWeightsPtr G3SkyMapWeights::Inv() const { G3SkyMapWeightsPtr out = this->Clone(false); out->TT->ConvertToDense(); - if (IsPolarized()) { + if (!!TQ) out->TQ->ConvertToDense(); + if (!!TU) out->TU->ConvertToDense(); + if (!!QQ) out->QQ->ConvertToDense(); + if (!!QU) out->QU->ConvertToDense(); + if (!!UU) out->UU->ConvertToDense(); - } for (size_t i = 0; i < TT->size(); i++) (*out)[i] = this->at(i).Inv(); @@ -856,20 +1300,12 @@ G3SkyMapWeightsPtr G3SkyMapWeights::Rebin(size_t scale) const g3_assert(IsCongruent()); G3SkyMapWeightsPtr out(new G3SkyMapWeights()); - out->TT = TT->Rebin(scale, false); - if (IsPolarized()) { - out->TQ = TQ->Rebin(scale, false); - out->TU = TU->Rebin(scale, false); - out->QQ = QQ->Rebin(scale, false); - out->QU = QU->Rebin(scale, false); - out->UU = UU->Rebin(scale, false); - } else { - out->TQ = NULL; - out->TU = NULL; - out->QQ = NULL; - out->QU = NULL; - out->UU = NULL; - } + out->TT = !TT ? NULL : TT->Rebin(scale, false); + out->TQ = !TQ ? NULL : TQ->Rebin(scale, false); + out->TU = !TU ? NULL : TU->Rebin(scale, false); + out->QQ = !QQ ? NULL : QQ->Rebin(scale, false); + out->QU = !QU ? NULL : QU->Rebin(scale, false); + out->UU = !UU ? NULL : UU->Rebin(scale, false); return out; } @@ -878,14 +1314,18 @@ void G3SkyMapWeights::Compact(bool zero_nans) { g3_assert(IsCongruent()); - TT->Compact(zero_nans); - if (IsPolarized()) { + if (!!TT) + TT->Compact(zero_nans); + if (!!TQ) TQ->Compact(zero_nans); + if (!!TU) TU->Compact(zero_nans); + if (!!QQ) QQ->Compact(zero_nans); + if (!!QU) QU->Compact(zero_nans); + if (!!UU) UU->Compact(zero_nans); - } } PYBINDINGS("maps") { @@ -960,6 +1400,7 @@ PYBINDINGS("maps") { .def("compatible", &G3SkyMap::IsCompatible, "Returns true if the input argument is a map with matching dimensions " "and boundaries on the sky.") + .def("nonzero", &pyskymap_nonzero, "Return indices of non-zero pixels in the map") .def("angles_to_pixels", &G3SkyMap::AnglesToPixels, (bp::arg("alphas"), bp::arg("deltas")), @@ -1038,6 +1479,18 @@ PYBINDINGS("maps") { "that is already sparse will be compactified in place in its " "current representation without additional memory overhead.") + .def("to_mask", &G3SkyMap::MakeMask, + (bp::arg("zero_nans")=false, bp::arg("zero_infs")=false), + "Create a G3SkyMapMask object from the parent map, with pixels " + "set to true where the map is non-zero (and optionally non-nan " + "and/or finite).") + + .def("apply_mask", &G3SkyMap::ApplyMask, + (bp::arg("mask"), bp::arg("inverse")=false), + "Apply a mask in-place to the map, optionally inverting which " + "pixels are zeroed. If inverse = False, this is equivalent to " + "in-place multiplication by the mask.") + .def(bp::self += bp::self) .def(bp::self *= bp::self) .def(bp::self -= bp::self) @@ -1053,8 +1506,11 @@ PYBINDINGS("maps") { .def("__sub__", &pyskymap_sub) .def("__sub__", &pyskymap_subd) .def("__rsub__", &pyskymap_rsubd) + .def("__imul__", &pyskymap_imultm) .def("__mul__", &pyskymap_mult) + .def("__mul__", &pyskymap_multm) .def("__mul__", &pyskymap_multd) + .def("__rmul__", &pyskymap_multm) .def("__rmul__", &pyskymap_multd) .def("__div__", &pyskymap_div) .def("__div__", &pyskymap_divd) @@ -1067,22 +1523,51 @@ PYBINDINGS("maps") { .def("__pow__", &pyskymap_pow) .def("__pow__", &pyskymap_powd) - .def("__lt__", &pyskymap_lt) - .def("__lt__", &pyskymap_ltd) - .def("__le__", &pyskymap_le) - .def("__le__", &pyskymap_led) - .def("__eq__", &pyskymap_eq) - .def("__eq__", &pyskymap_eqd) - .def("__ne__", &pyskymap_ne) - .def("__ne__", &pyskymap_ned) - .def("__ge__", &pyskymap_ge) - .def("__ge__", &pyskymap_ged) - .def("__gt__", &pyskymap_gt) - .def("__gt__", &pyskymap_gtd) + .def(bp::self < bp::self) + .def(bp::self <= bp::self) + .def(bp::self == bp::self) + .def(bp::self != bp::self) + .def(bp::self >= bp::self) + .def(bp::self > bp::self) + .def(bp::self < double()) + .def(bp::self <= double()) + .def(bp::self == double()) + .def(bp::self != double()) + .def(bp::self >= double()) + .def(bp::self > double()) + .def("__bool__", &pyskymap_bool) - .def("any", &pyskymap_any) - .def("all", &pyskymap_all) + .def("_cany", &G3SkyMap::any, (bp::arg("where")=G3SkyMapMaskConstPtr())) + .def("_call", &G3SkyMap::all, (bp::arg("where")=G3SkyMapMaskConstPtr())) + .def("_csum", &G3SkyMap::sum, (bp::arg("where")=G3SkyMapMaskConstPtr())) + .def("_cmean", &G3SkyMap::mean, (bp::arg("where")=G3SkyMapMaskConstPtr())) + .def("median", &G3SkyMap::median, (bp::arg("where")=G3SkyMapMaskConstPtr())) + .def("_cvar", &G3SkyMap::var, (bp::arg("ddof")=0, + bp::arg("where")=G3SkyMapMaskConstPtr())) + .def("_cstd", &G3SkyMap::std, (bp::arg("ddof")=0, + bp::arg("where")=G3SkyMapMaskConstPtr())) + .def("_cmin", &G3SkyMap::min, (bp::arg("where")=G3SkyMapMaskConstPtr())) + .def("_cmax", &G3SkyMap::max, (bp::arg("where")=G3SkyMapMaskConstPtr())) + .def("_cargmin", &G3SkyMap::argmin, (bp::arg("where")=G3SkyMapMaskConstPtr())) + .def("_cargmax", &G3SkyMap::argmax, (bp::arg("where")=G3SkyMapMaskConstPtr())) + + .def("nansum", &G3SkyMap::nansum, (bp::arg("where")=G3SkyMapMaskConstPtr())) + .def("nanmean", &G3SkyMap::nanmean, (bp::arg("where")=G3SkyMapMaskConstPtr())) + .def("nanmedian", &G3SkyMap::nanmedian, (bp::arg("where")=G3SkyMapMaskConstPtr())) + .def("nanvar", &G3SkyMap::nanvar, (bp::arg("ddof")=0, + bp::arg("where")=G3SkyMapMaskConstPtr())) + .def("nanstd", &G3SkyMap::nanstd, (bp::arg("ddof")=0, + bp::arg("where")=G3SkyMapMaskConstPtr())) + .def("nanmin", &G3SkyMap::nanmin, (bp::arg("where")=G3SkyMapMaskConstPtr())) + .def("nanmax", &G3SkyMap::nanmax, (bp::arg("where")=G3SkyMapMaskConstPtr())) + .def("nanargmin", &G3SkyMap::nanargmin, (bp::arg("where")=G3SkyMapMaskConstPtr())) + .def("nanargmax", &G3SkyMap::nanargmax, (bp::arg("where")=G3SkyMapMaskConstPtr())) + + .def("isinf", &G3SkyMap::isinf, (bp::arg("where")=G3SkyMapMaskConstPtr())) + .def("isnan", &G3SkyMap::isnan, (bp::arg("where")=G3SkyMapMaskConstPtr())) + .def("isfinite", &G3SkyMap::isfinite, (bp::arg("where")=G3SkyMapMaskConstPtr())) ; + boost::python::implicitly_convertible(); EXPORT_FRAMEOBJECT(G3SkyMapWeights, init<>(), "Polarized (Mueller matrix) or unpolarized (scalar) map pixel weights.") @@ -1120,6 +1605,12 @@ PYBINDINGS("maps") { .def("inv", &G3SkyMapWeights::Inv, "Return the inverse of the Mueller matrix for each pixel") + .def("apply_mask", &G3SkyMapWeights::ApplyMask, + (bp::arg("mask"), bp::arg("inverse")=false), + "Apply a mask in-place to the weights, optionally inverting which " + "pixels are zeroed. If inverse = False, this is equivalent to " + "in-place multiplication by the mask.") + .def("clone", &G3SkyMapWeights::Clone, (bp::arg("copy_data")=true), "Return weights of the same type, populated with a copy of the data " "if the argument is true (default), empty otherwise.") @@ -1131,8 +1622,12 @@ PYBINDINGS("maps") { .def(bp::self /= double()) .def("__add__", &pyskymapweights_add) + .def("__imul__", &pyskymapweights_imultma) .def("__mul__", &pyskymapweights_multm) + .def("__mul__", &pyskymapweights_multma) .def("__mul__", &pyskymapweights_multd) + .def("__rmul__", &pyskymapweights_multm) + .def("__rmul__", &pyskymapweights_multma) .def("__rmul__", &pyskymapweights_multd) .def("__div__", &pyskymapweights_divd) .def("__truediv__", &pyskymapweights_divd) diff --git a/maps/src/G3SkyMapMask.cxx b/maps/src/G3SkyMapMask.cxx new file mode 100644 index 00000000..ece15b9e --- /dev/null +++ b/maps/src/G3SkyMapMask.cxx @@ -0,0 +1,594 @@ +#include +#include +#include +#include + +#include +#include +#include + +G3SkyMapMask::G3SkyMapMask(const G3SkyMap &parent, bool use_data, + bool zero_nans, bool zero_infs) : G3FrameObject() +{ + // masks have no units + G3SkyMapPtr tmp = parent.Clone(false); + tmp->units = G3Timestream::None; + tmp->pol_type = G3SkyMap::None; + tmp->SetPolConv(G3SkyMap::ConvNone); + tmp->weighted = false; + parent_ = tmp; + data_ = std::vector(parent.size()); + + if (use_data) + FillFromMap(parent, zero_nans, zero_infs); +} + +G3SkyMapMask::G3SkyMapMask(const G3SkyMap &parent, boost::python::object v, + bool zero_nans, bool zero_infs) : G3FrameObject() +{ + // masks have no units + G3SkyMapPtr tmp = parent.Clone(false); + tmp->units = G3Timestream::None; + tmp->pol_type = G3SkyMap::None; + tmp->SetPolConv(G3SkyMap::ConvNone); + tmp->weighted = false; + parent_ = tmp; + data_ = std::vector(parent.size()); + + // fall back to simple constructor + if (boost::python::extract(v).check()) { + if (boost::python::extract(v)()) + FillFromMap(parent, zero_nans, zero_infs); + return; + } + + FillFromArray(v, zero_nans, zero_infs); +} + +G3SkyMapMask::G3SkyMapMask(const G3SkyMapMask &m) : G3FrameObject() +{ + parent_ = m.parent_; + data_ = std::vector(m.data_); +} + +G3SkyMapMaskPtr +G3SkyMapMask::Clone(bool copy_data) const +{ + if (copy_data) + return boost::make_shared(*this); + else + return boost::make_shared(*Parent()); +} + +void +G3SkyMapMask::FillFromMap(const G3SkyMap &m, bool zero_nans, bool zero_infs) +{ + g3_assert(IsCompatible(m)); + + for (size_t i = 0; i < m.size(); i++) { + double v = m.at(i); + if (v == 0) + continue; + if (zero_nans && std::isnan(v)) + continue; + if (zero_infs && std::isinf(v)) + continue; + data_[i] = true; + } +} + +#define FILL_BUFFER(type) \ + for (size_t i = 0; i < view.len / sizeof(type); i++) { \ + double v = ((type *)view.buf)[i]; \ + if (v == 0) \ + continue; \ + if (zero_nans && std::isnan(v)) \ + continue; \ + if (zero_infs && std::isinf(v)) \ + continue; \ + data_[i] = true; \ + } + +void +G3SkyMapMask::FillFromArray(boost::python::object v, bool zero_nans, bool zero_infs) +{ + + Py_buffer view; + + if (PyObject_GetBuffer(v.ptr(), &view, + PyBUF_FORMAT | PyBUF_C_CONTIGUOUS) != -1) { + if (view.ndim != 1) { + PyBuffer_Release(&view); + log_fatal("Only 1-D masks supported"); + } + + size_t npix = view.shape[0]; + if (npix != size()) { + PyBuffer_Release(&view); + log_fatal("Got array of shape (%zu,), expected (%zu,)", npix, size()); + } + + const char *format = view.format; + + if (format[0] == '@' || format[0] == '=') + format++; +#if BYTE_ORDER == LITTLE_ENDIAN + else if (format[0] == '<') + format++; + else if (format[0] == '>' || format[0] == '!') { + PyBuffer_Release(&view); + log_fatal("Does not support big-endian numpy arrays"); + } +#else + else if (format[0] == '<') { + PyBuffer_Release(&view); + log_fatal("Does not support little-endian numpy arrays"); + } else if (format[0] == '>' || format[0] == '!') + format++; +#endif + + if (strcmp(format, "d") == 0) { + FILL_BUFFER(double); + } else if (strcmp(format, "f") == 0) { + FILL_BUFFER(float); + } else if (strcmp(format, "i") == 0) { + FILL_BUFFER(int); + } else if (strcmp(format, "I") == 0) { + FILL_BUFFER(unsigned int); + } else if (strcmp(format, "l") == 0) { + FILL_BUFFER(long); + } else if (strcmp(format, "L") == 0) { + FILL_BUFFER(unsigned long); + } else if (strcmp(format, "b") == 0) { + FILL_BUFFER(char); + } else if (strcmp(format, "B") == 0) { + FILL_BUFFER(unsigned char); + } else if (strcmp(format, "?") == 0) { + FILL_BUFFER(bool); + } else { + PyBuffer_Release(&view); + log_fatal("Unknown type code %s", view.format); + } + PyBuffer_Release(&view); + + return; + } + + throw boost::python::error_already_set(); +} + +G3SkyMapMaskPtr +G3SkyMapMask::ArrayClone(boost::python::object v, bool zero_nans, bool zero_infs) const +{ + G3SkyMapMaskPtr m = Clone(false); + m->FillFromArray(v, zero_nans, zero_infs); + return m; +} + +G3SkyMapMask::const_iterator::const_iterator(const G3SkyMapMask &mask, bool begin) : + mask_(mask) +{ + index_ = begin ? 0 : mask_.size(); + set_value(); +} + +G3SkyMapMask::const_iterator +G3SkyMapMask::const_iterator::operator++() +{ + ++index_; + set_value(); + + return *this; +} + +size_t +G3SkyMapMask::size() const +{ + return data_.size(); +} + +std::vector::reference +G3SkyMapMask::operator [] (size_t i) +{ + return data_[i]; +} + +bool +G3SkyMapMask::at (size_t i) const +{ + if (i >= data_.size()) + return false; + return data_[i]; +} + +G3SkyMapMask & +G3SkyMapMask::operator |=(const G3SkyMapMask &rhs) +{ + g3_assert(IsCompatible(rhs)); + + for (size_t i = 0; i < size(); i++) + (*this)[i] = rhs.at(i) || at(i); + + return *this; +} + +G3SkyMapMask & +G3SkyMapMask::operator &=(const G3SkyMapMask &rhs) +{ + g3_assert(IsCompatible(rhs)); + + for (auto i: *this) + (*this)[i.first] = rhs.at(i.first) && i.second; + + return *this; +} + +G3SkyMapMask & +G3SkyMapMask::operator ^=(const G3SkyMapMask &rhs) +{ + g3_assert(IsCompatible(rhs)); + + for (size_t i = 0; i < size(); i++) + (*this)[i] = rhs.at(i) ^ at(i); + + return *this; +} + +G3SkyMapMask & +G3SkyMapMask::invert() +{ + for (size_t i = 0; i < size(); i++) + (*this)[i] = !at(i); + + return *this; +} + +bool +G3SkyMapMask::all() const +{ + for (size_t i = 0; i < size(); i++) + if (at(i) == 0) + return false; + return true; +} + +bool +G3SkyMapMask::any() const +{ + for (auto i: *this) + if (i.second != 0) + return true; + return false; +} + +size_t +G3SkyMapMask::sum() const +{ + size_t s = 0; + for (auto i: *this) + s += i.second; + return s; +} + +std::vector +G3SkyMapMask::NonZeroPixels() const +{ + std::vector indices; + + for (auto i: *this) { + if (i.second) + indices.push_back(i.first); + } + + return indices; +} + +G3SkyMapMask +G3SkyMapMask::operator ~() const +{ + G3SkyMapMask mask(*Parent()); + for (size_t i = 0; i < size(); i++) + if (!at(i)) + mask[i] = true; + + return mask; +} + +G3SkyMapMask +G3SkyMapMask::operator |(const G3SkyMapMask &rhs) const +{ + g3_assert(IsCompatible(rhs)); + + G3SkyMapMask mask(*Parent()); + for (size_t i = 0; i < size(); i++) + if (at(i) || rhs.at(i)) + mask[i] = true; + + return mask; +} + +G3SkyMapMask +G3SkyMapMask::operator &(const G3SkyMapMask &rhs) const +{ + g3_assert(IsCompatible(rhs)); + + G3SkyMapMask mask(*Parent()); + for (auto i: *this) + if (i.second && rhs.at(i.first)) + mask[i.first] = true; + + return mask; +} + +G3SkyMapMask +G3SkyMapMask::operator ^(const G3SkyMapMask &rhs) const +{ + g3_assert(IsCompatible(rhs)); + + G3SkyMapMask mask(*Parent()); + for (size_t i = 0; i < size(); i++) + if (at(i) ^ rhs.at(i)) + mask[i] = true; + + return mask; +} + +G3SkyMapMask +G3SkyMapMask::operator ==(const G3SkyMapMask &rhs) const +{ + g3_assert(IsCompatible(rhs)); + + G3SkyMapMask mask(*Parent()); + for (size_t i = 0; i < size(); i++) + if (at(i) == rhs.at(i)) + mask[i] = true; + + return mask; +} + +G3SkyMapMask +G3SkyMapMask::operator !=(const G3SkyMapMask &rhs) const +{ + g3_assert(IsCompatible(rhs)); + + G3SkyMapMask mask(*Parent()); + for (size_t i = 0; i < size(); i++) + if (at(i) != rhs.at(i)) + mask[i] = true; + + return mask; +} + +void +G3SkyMapMask::ApplyMask(const G3SkyMapMask &rhs, bool inverse) +{ + g3_assert(IsCompatible(rhs)); + + for (auto i: *this) + if (i.second && rhs.at(i.first) == inverse) + (*this)[i.first] = false; +} + +bool +G3SkyMapMask::IsCompatible(const G3SkyMap &map) const +{ + return Parent()->IsCompatible(map); +} + +bool +G3SkyMapMask::IsCompatible(const G3SkyMapMask &mask) const +{ + return Parent()->IsCompatible(*mask.Parent()); +} + +G3SkyMapPtr +G3SkyMapMask::MakeBinaryMap() const +{ + G3SkyMapPtr out = Parent()->Clone(); + + for (auto i: *this) { + if (i.second) + (*out)[i.first] = 1.0; + } + + return out; +} + +template +void G3SkyMapMask::serialize(A &ar, unsigned v) +{ + using namespace cereal; + ar & make_nvp("G3FrameObject", base_class(this)); + ar & make_nvp("parent", parent_); + ar & make_nvp("data", data_); +} + +G3_SERIALIZABLE_CODE(G3SkyMapMask); + +static bool +skymapmask_getitem(const G3SkyMapMask &m, boost::python::object index) +{ + using namespace boost::python; + + int i = 0; + if (extract(index).check()) { + i = extract(index)(); + + if (i < 0) + i = m.size() + i; + } else if (extract(index).check()) { + int x, y; + + tuple t = extract(index)(); + FlatSkyMapConstPtr fsm = boost::dynamic_pointer_cast(m.Parent()); + if (!fsm) { + PyErr_SetString(PyExc_TypeError, + "N-D pixels, but underlying map is not a flat sky map"); + boost::python::throw_error_already_set(); + } + + x = extract(t[1])(); + y = extract(t[0])(); + if (x < 0) + x += fsm->shape()[0]; + if (y < 0) + y += fsm->shape()[0]; + if (size_t(x) >= fsm->shape()[0] || + size_t(y) >= fsm->shape()[1]) { + PyErr_SetString(PyExc_IndexError, "Index out of range"); + boost::python::throw_error_already_set(); + } + + i = y * fsm->shape()[0] + x; + } else { + PyErr_SetString(PyExc_TypeError, + "Need to pass an integer pixel ID or (optionally) for 2D maps a tuple of coordinates"); + boost::python::throw_error_already_set(); + } + + if (i < 0 || size_t(i) >= m.size()) { + PyErr_SetString(PyExc_IndexError, "Index out of range"); + boost::python::throw_error_already_set(); + } + + return m.at(i); +} + +static void +skymapmask_setitem(G3SkyMapMask &m, boost::python::object index, bool val) +{ + using namespace boost::python; + + int i = 0; + if (extract(index).check()) { + i = extract(index)(); + + if (i < 0) + i = m.size() + i; + } else if (extract(index).check()) { + int x, y; + + tuple t = extract(index)(); + FlatSkyMapConstPtr fsm = boost::dynamic_pointer_cast(m.Parent()); + if (!fsm) { + PyErr_SetString(PyExc_TypeError, + "N-D pixels, but underlying map is not a flat sky map"); + boost::python::throw_error_already_set(); + } + + x = extract(t[1])(); + y = extract(t[0])(); + if (x < 0) + x += fsm->shape()[0]; + if (y < 0) + y += fsm->shape()[0]; + if (size_t(x) >= fsm->shape()[0] || + size_t(y) >= fsm->shape()[1]) { + PyErr_SetString(PyExc_IndexError, "Index out of range"); + boost::python::throw_error_already_set(); + } + + i = y * fsm->shape()[0] + x; + } else { + PyErr_SetString(PyExc_TypeError, + "Need to pass an integer pixel ID or (optionally) for 2D maps a tuple of coordinates"); + boost::python::throw_error_already_set(); + } + + if (i < 0 || size_t(i) >= m.size()) { + PyErr_SetString(PyExc_IndexError, "Index out of range"); + boost::python::throw_error_already_set(); + } + + m[i] = val; +} + +static G3SkyMapMaskPtr +skymapmask_pyinvert(G3SkyMapMaskPtr m) +{ + // Reference-counting problems returning references to Python, + // so use shared pointers. + m->invert(); + return m; +} + +// This function handles some implicit pointer conversions that boost won't +// do because G3SkyMap is an abstract type. Just do it by hand. +static G3SkyMapPtr +skymapmask_pyparent(G3SkyMapMaskPtr m) +{ + return boost::const_pointer_cast(m->Parent()); +} + +static bool +skymapmask_pybool(G3SkyMapMaskPtr m) +{ + PyErr_SetString(PyExc_ValueError, + "ValueError: The truth value of a G3SkyMapMask is ambiguous. Use m.any() or m.all()"); + bp::throw_error_already_set(); + + return false; +} + +PYBINDINGS("maps") +{ + using namespace boost::python; + + EXPORT_FRAMEOBJECT(G3SkyMapMask, no_init, + "Boolean mask of a sky map. Set pixels to use to true, pixels to " + "ignore to false. If use_data set in contrast, mask initialized to " + "true where input map is non-zero; otherwise, all elements are " + "initialized to zero. Use zero_nans or zero_infs to exclude nan " + "or inf elements from the mask.") + .def(bp::init( + (bp::arg("parent"), + bp::arg("use_data")=false, + bp::arg("zero_nans")=false, + bp::arg("zero_infs")=false), + "Instantiate a G3SkyMapMask from a parent G3SkyMap")) + .def(bp::init( + (bp::arg("parent"), + bp::arg("data"), + bp::arg("zero_nans")=false, + bp::arg("zero_infs")=false), + "Instantiate a G3SkyMapMask from a 1D numpy array")) + .def("clone", &G3SkyMapMask::Clone, (bp::arg("copy_data")=true), + "Return a mask of the same type, populated with a copy of the data " + "if the argument is true (default), empty otherwise.") + .def("array_clone", &G3SkyMapMask::ArrayClone, + (bp::arg("data"), bp::arg("zero_nans")=false, bp::arg("zero_infs")=false), + "Return a mask of the same type, populated from the input numpy array") + .add_property("parent", &skymapmask_pyparent, "\"Parent\" map which " + "contains no data, but can be used to retrieve the parameters of " + "the map to which this mask corresponds.") + .add_property("size", &G3SkyMapMask::size, "Number of pixels in mask") + .def("compatible", &G3SkyMapMask::IsCompatible, "Returns true if the two masks can be applied to the same map.") + .def("compatible", &G3SkyMapMask::IsCompatible, "Returns true if this mask can be applied to the given map.") + .def("__getitem__", &skymapmask_getitem) + .def("__setitem__", &skymapmask_setitem) + .def("__bool__", &skymapmask_pybool) + .def("invert", &skymapmask_pyinvert, "Invert all elements in mask") + .def("_call", &G3SkyMapMask::all, "Test whether all elements are non-zero") + .def("_cany", &G3SkyMapMask::any, "Test whether any elements are non-zero") + .def("_csum", &G3SkyMapMask::sum, "Sum of all elements in mask") + .def("nonzero", &G3SkyMapMask::NonZeroPixels, + "Return a list of indices of non-zero pixels in the mask") + .def("apply_mask", &G3SkyMapMask::ApplyMask, + (bp::arg("mask"), bp::arg("inverse")=false), + "Apply a mask in-place to the mask, optionally inverting which " + "pixels are zeroed. If inverse = False, this is equivalent to " + "in-place element-wise logical-and with the mask.") + .def(bp::self |= bp::self) + .def(bp::self &= bp::self) + .def(bp::self ^= bp::self) + .def(~bp::self) + .def(bp::self | bp::self) + .def(bp::self & bp::self) + .def(bp::self ^ bp::self) + .def(bp::self == bp::self) + .def(bp::self != bp::self) + .def("to_map", &G3SkyMapMask::MakeBinaryMap, "Create a skymap with data set to the contents of this mask (1.0 where True, 0.0 where False), which can be useful for plotting.") + ; + register_pointer_conversions(); +} + diff --git a/maps/src/HealpixSkyMap.cxx b/maps/src/HealpixSkyMap.cxx index b08dafb2..3d006770 100644 --- a/maps/src/HealpixSkyMap.cxx +++ b/maps/src/HealpixSkyMap.cxx @@ -9,6 +9,7 @@ #include #include #include +#include #include "mapdata.h" @@ -104,7 +105,7 @@ HealpixSkyMap::HealpixSkyMap(const HealpixSkyMap & fm) : if (fm.dense_) dense_ = new std::vector(*fm.dense_); else if (fm.ring_sparse_) - ring_sparse_ = new SparseMapData(*fm.ring_sparse_); + ring_sparse_ = new SparseMapData(*fm.ring_sparse_); else if (fm.indexed_sparse_) indexed_sparse_ = new std::unordered_map( *fm.indexed_sparse_); @@ -182,7 +183,7 @@ HealpixSkyMap::load(A &ar, unsigned v) ar & make_nvp("data", *dense_); break; case 2: - ring_sparse_ = new SparseMapData(1,1); + ring_sparse_ = new SparseMapData(1,1); ar & make_nvp("data", *ring_sparse_); break; case 1: @@ -263,7 +264,7 @@ HealpixSkyMap::const_iterator::operator++() ++index_; ++it_dense_; } else if (map_.ring_sparse_) { - SparseMapData::const_iterator iter(*map_.ring_sparse_, j_, k_); + SparseMapData::const_iterator iter(*map_.ring_sparse_, j_, k_); ++iter; j_ = iter.x; k_ = iter.y; @@ -308,7 +309,7 @@ HealpixSkyMap::ConvertToRingSparse() if (ring_sparse_) return; - ring_sparse_ = new SparseMapData(info_.nring(), info_.nring()); + ring_sparse_ = new SparseMapData(info_.nring(), info_.nring()); if (dense_) { auto *old_dense = dense_; @@ -614,7 +615,7 @@ HealpixSkyMap::operator [] (size_t i) return (*indexed_sparse_)[i]; if (!ring_sparse_) - ring_sparse_ = new SparseMapData(info_.nring(), info_.nring()); + ring_sparse_ = new SparseMapData(info_.nring(), info_.nring()); auto ridx = info_.PixelToRing(i); return (*ring_sparse_)(ridx.first, ridx.second); @@ -698,6 +699,19 @@ G3SkyMap &HealpixSkyMap::operator *=(const G3SkyMap &rhs) { return *this; } +G3SkyMap & +HealpixSkyMap::operator *=(const G3SkyMapMask &rhs) +{ + g3_assert(rhs.IsCompatible(*this)); + + for (auto i : *this) { + if (!rhs.at(i.first) && i.second != 0) + (*this)[i.first] = 0; + } + + return *this; +} + G3SkyMap &HealpixSkyMap::operator /=(const G3SkyMap &rhs) { g3_assert(IsCompatible(rhs)); if (units == G3Timestream::None) @@ -930,6 +944,16 @@ HealpixSkyMap::NonZeroPixels(std::vector &indices, } } +void +HealpixSkyMap::ApplyMask(const G3SkyMapMask &mask, bool inverse) +{ + g3_assert(mask.IsCompatible(*this)); + + for (auto i: *this) + if (i.second != 0 && mask.at(i.first) == inverse) + (*this)[i.first] = 0; +} + std::vector HealpixSkyMap::shape() const @@ -1055,7 +1079,7 @@ void HealpixSkyMap::SetShiftRa(bool shift) HealpixSkyMapInfo info(info_); info.SetShifted(shift); - SparseMapData *ring_sparse = new SparseMapData(info_.nring(), info_.nring()); + SparseMapData *ring_sparse = new SparseMapData(info_.nring(), info_.nring()); for (auto i : *this) { if (i.second != 0) { auto ridx = info.PixelToRing(i.first); @@ -1124,6 +1148,73 @@ HealpixSkyMap_nonzeropixels(const HealpixSkyMap &m) return boost::python::make_tuple(i, d); } +static double +skymap_getitem(const G3SkyMap &skymap, int i) +{ + + if (i < 0) + i = skymap.size() + i; + if (size_t(i) >= skymap.size()) { + PyErr_SetString(PyExc_IndexError, "Index out of range"); + bp::throw_error_already_set(); + } + + return skymap.at(i); +} + +static void +skymap_setitem(G3SkyMap &skymap, int i, double val) +{ + + if (i < 0) + i = skymap.size() + i; + if (size_t(i) >= skymap.size()) { + PyErr_SetString(PyExc_IndexError, "Index out of range"); + bp::throw_error_already_set(); + } + + skymap[i] = val; +} + +static std::vector +HealpixSkyMap_getitem_masked(const HealpixSkyMap &skymap, const G3SkyMapMask &m) +{ + g3_assert(m.IsCompatible(skymap)); + std::vector out; + + for (auto i : skymap) { + if (m.at(i.first)) + out.push_back(i.second); + } + + return out; +} + +static void +HealpixSkyMap_setitem_masked(HealpixSkyMap &skymap, const G3SkyMapMask &m, + bp::object val) +{ + g3_assert(m.IsCompatible(skymap)); + + if (bp::extract(val).check()) { + double dval = bp::extract(val)(); + for (auto i : skymap) { + if (m.at(i.first)) + skymap[i.first] = dval; + } + } else { + // XXX: the iterable case probably be optimized for numpy arrays + // XXX: check for size congruence first? + size_t j = 0; + for (auto i : skymap) { + if (m.at(i.first)) { + skymap[i.first] = bp::extract(val[j])(); + j++; + } + } + } +} + static int HealpixSkyMap_getbuffer(PyObject *obj, Py_buffer *view, int flags) { @@ -1266,6 +1357,11 @@ PYBINDINGS("maps") "holes or very small filling factors. " "If set to True, converts the map to this representation." ) + .def("__getitem__", &skymap_getitem) + .def("__setitem__", &skymap_setitem) + .def("__getitem__", HealpixSkyMap_getitem_masked) + .def("__setitem__", HealpixSkyMap_setitem_masked) + .def("nonzero_pixels", &HealpixSkyMap_nonzeropixels, "Returns a list of the indices of the non-zero pixels in the " "map and a list of the values of those non-zero pixels.") @@ -1284,8 +1380,6 @@ PYBINDINGS("maps") #endif implicitly_convertible(); - implicitly_convertible(); implicitly_convertible(); - implicitly_convertible(); } diff --git a/maps/src/MapTODMasker.cxx b/maps/src/MapTODMasker.cxx index 8fbc0e6f..0c956ca6 100644 --- a/maps/src/MapTODMasker.cxx +++ b/maps/src/MapTODMasker.cxx @@ -15,7 +15,7 @@ class MapTODMasker : public G3Module { public: MapTODMasker(std::string pointing, std::string timestreams, - G3SkyMapConstPtr mask, std::string tod_mask, + G3SkyMapMaskConstPtr mask, std::string tod_mask, std::string bolo_properties_name); virtual ~MapTODMasker() {} @@ -24,7 +24,7 @@ class MapTODMasker : public G3Module { private: std::string pointing_; std::string timestreams_; - G3SkyMapConstPtr mask_; + G3SkyMapMaskConstPtr mask_; std::string output_; std::string boloprops_name_; @@ -34,7 +34,7 @@ class MapTODMasker : public G3Module { }; EXPORT_G3MODULE("maps", MapTODMasker, - (init((arg("pointing"), arg("timestreams"), arg("mask"), arg("tod_mask")="FilterMask", arg("bolo_properties_name")="BolometerProperties"))), @@ -50,7 +50,7 @@ EXPORT_G3MODULE("maps", MapTODMasker, "can generally be left at its default value."); MapTODMasker::MapTODMasker(std::string pointing, std::string timestreams, - G3SkyMapConstPtr mask, std::string tod_mask, + G3SkyMapMaskConstPtr mask, std::string tod_mask, std::string bolo_properties_name) : pointing_(pointing), timestreams_(timestreams), mask_(mask), output_(tod_mask), boloprops_name_(bolo_properties_name) @@ -81,6 +81,8 @@ MapTODMasker::Process(G3FramePtr frame, std::deque &out) return; } + G3SkyMapConstPtr skymap = mask_->Parent(); + G3TimestreamMapConstPtr tsm = frame->Get(timestreams_); G3MapVectorBoolPtr output(new G3MapVectorBool); @@ -105,12 +107,12 @@ MapTODMasker::Process(G3FramePtr frame, std::deque &out) // Get per-detector pointing timestream std::vector alpha, delta; get_detector_pointing(bp.x_offset, bp.y_offset, *pointing, - mask_->coord_ref, alpha, delta); - auto detpointing = mask_->AnglesToPixels(alpha, delta); + skymap->coord_ref, alpha, delta); + auto detpointing = skymap->AnglesToPixels(alpha, delta); det.resize(pointing->size()); for (size_t j = 0; j < det.size(); j++) - det[j] = !!(mask_->at(detpointing[j])); + det[j] = mask_->at(detpointing[j]); // Find out if any elements are set, delete entry if not bool is_set = false; diff --git a/maps/src/mapdata.cxx b/maps/src/mapdata.cxx index f824774b..ac2dac2a 100644 --- a/maps/src/mapdata.cxx +++ b/maps/src/mapdata.cxx @@ -3,8 +3,9 @@ #include "mapdata.h" -SparseMapData::const_iterator -SparseMapData::const_iterator::operator++() { +template +typename SparseMapData::const_iterator +SparseMapData::const_iterator::operator++() { const_iterator end(sparse_.end()); if (x > end.x || sparse_.data_.size() == 0) { @@ -19,7 +20,7 @@ SparseMapData::const_iterator::operator++() { return *this; } - const SparseMapData::data_element &column = sparse_.data_[x - sparse_.offset_]; + const SparseMapData::data_element &column = sparse_.data_[x - sparse_.offset_]; if (column.second.size() > 0 && y < column.first) { y = column.first; return *this; @@ -43,8 +44,9 @@ SparseMapData::const_iterator::operator++() { return *this; } +template DenseMapData * -SparseMapData::to_dense() const +SparseMapData::to_dense() const { DenseMapData *rv = new DenseMapData(xlen_, ylen_); for (size_t ix = 0; ix < data_.size(); ix++) { @@ -58,7 +60,8 @@ SparseMapData::to_dense() const return rv; } -SparseMapData::SparseMapData(const DenseMapData &dense_map) : +template +SparseMapData::SparseMapData(const DenseMapData &dense_map) : xlen_(dense_map.xdim()), ylen_(dense_map.ydim()), offset_(0) { double val; @@ -72,8 +75,9 @@ SparseMapData::SparseMapData(const DenseMapData &dense_map) : } } +template void -SparseMapData::compact() +SparseMapData::compact() { if (data_.size() == 0) return; @@ -104,8 +108,9 @@ SparseMapData::compact() } -SparseMapData & -SparseMapData::operator+=(const SparseMapData &r) +template <> +SparseMapData & +SparseMapData::operator+=(const SparseMapData &r) { assert(xlen_ == r.xlen_); assert(ylen_ == r.ylen_); @@ -124,8 +129,9 @@ SparseMapData::operator+=(const SparseMapData &r) return *this; } -SparseMapData & -SparseMapData::operator+=(const DenseMapData &r) +template <> +SparseMapData & +SparseMapData::operator+=(const DenseMapData &r) { assert(xlen_ == r.xdim()); assert(ylen_ == r.ydim()); @@ -141,8 +147,9 @@ SparseMapData::operator+=(const DenseMapData &r) return *this; } -SparseMapData & -SparseMapData::operator-=(const SparseMapData &r) +template <> +SparseMapData & +SparseMapData::operator-=(const SparseMapData &r) { assert(xlen_ == r.xlen_); assert(ylen_ == r.ylen_); @@ -161,8 +168,9 @@ SparseMapData::operator-=(const SparseMapData &r) return *this; } -SparseMapData & -SparseMapData::operator-=(const DenseMapData &r) +template <> +SparseMapData & +SparseMapData::operator-=(const DenseMapData &r) { assert(xlen_ == r.xdim()); assert(ylen_ == r.ydim()); @@ -178,8 +186,9 @@ SparseMapData::operator-=(const DenseMapData &r) return *this; } -SparseMapData & -SparseMapData::operator*=(const SparseMapData &r) +template <> +SparseMapData & +SparseMapData::operator*=(const SparseMapData &r) { assert(xlen_ == r.xlen_); assert(ylen_ == r.ylen_); @@ -196,8 +205,9 @@ SparseMapData::operator*=(const SparseMapData &r) return *this; } -SparseMapData & -SparseMapData::operator*=(const DenseMapData &r) +template <> +SparseMapData & +SparseMapData::operator*=(const DenseMapData &r) { assert(xlen_ == r.xdim()); assert(ylen_ == r.ydim()); @@ -214,8 +224,9 @@ SparseMapData::operator*=(const DenseMapData &r) return *this; } -SparseMapData & -SparseMapData::operator*=(double r) +template <> +SparseMapData & +SparseMapData::operator*=(double r) { for (size_t ix = 0; ix < data_.size(); ix++) { size_t x = offset_ + ix; @@ -229,8 +240,9 @@ SparseMapData::operator*=(double r) return *this; } -SparseMapData & -SparseMapData::operator/=(const SparseMapData &r) +template <> +SparseMapData & +SparseMapData::operator/=(const SparseMapData &r) { assert(xlen_ == r.xlen_); assert(ylen_ == r.ylen_); @@ -249,8 +261,9 @@ SparseMapData::operator/=(const SparseMapData &r) return *this; } -SparseMapData & -SparseMapData::operator/=(const DenseMapData &r) +template <> +SparseMapData & +SparseMapData::operator/=(const DenseMapData &r) { assert(xlen_ == r.xdim()); assert(ylen_ == r.ydim()); @@ -269,8 +282,9 @@ SparseMapData::operator/=(const DenseMapData &r) return *this; } -SparseMapData & -SparseMapData::operator/=(double r) +template <> +SparseMapData & +SparseMapData::operator/=(double r) { assert(r != 0); @@ -302,7 +316,7 @@ DenseMapData::operator+=(const DenseMapData &r) } DenseMapData & -DenseMapData::operator+=(const SparseMapData &r) +DenseMapData::operator+=(const SparseMapData &r) { assert(xlen_ == r.xdim()); assert(ylen_ == r.ydim()); @@ -347,7 +361,7 @@ DenseMapData::operator-=(const DenseMapData &r) } DenseMapData & -DenseMapData::operator-=(const SparseMapData &r) +DenseMapData::operator-=(const SparseMapData &r) { assert(xlen_ == r.xdim()); assert(ylen_ == r.ydim()); @@ -392,7 +406,7 @@ DenseMapData::operator*=(const DenseMapData &r) } DenseMapData & -DenseMapData::operator*=(const SparseMapData &r) +DenseMapData::operator*=(const SparseMapData &r) { assert(xlen_ == r.xdim()); assert(ylen_ == r.ydim()); @@ -434,7 +448,7 @@ DenseMapData::operator/=(const DenseMapData &r) } DenseMapData & -DenseMapData::operator/=(const SparseMapData &r) +DenseMapData::operator/=(const SparseMapData &r) { assert(xlen_ == r.xdim()); assert(ylen_ == r.ydim()); @@ -460,3 +474,6 @@ DenseMapData::operator/=(double r) return *this; } +template class SparseMapData; +template class SparseMapData; + diff --git a/maps/src/mapdata.h b/maps/src/mapdata.h index b620083c..1be172ff 100644 --- a/maps/src/mapdata.h +++ b/maps/src/mapdata.h @@ -3,10 +3,9 @@ #include #include -class SparseMapData; class DenseMapData; - +template class SparseMapData { public: SparseMapData(size_t xlen, size_t ylen) : @@ -34,7 +33,7 @@ class SparseMapData { return nz; } - double at(size_t x, size_t y) const { + T at(size_t x, size_t y) const { if (x < offset_ || x >= offset_ + data_.size()) return 0; const data_element &column = data_[x-offset_]; @@ -44,9 +43,9 @@ class SparseMapData { return column.second[y-column.first]; } - double operator()(size_t x, size_t y) const { return at(x, y); } + T operator()(size_t x, size_t y) const { return at(x, y); } - double &operator()(size_t x, size_t y) { + typename std::vector::reference operator()(size_t x, size_t y) { assert(x >= 0); assert(x < xlen_); assert(y >= 0); @@ -65,31 +64,31 @@ class SparseMapData { if (column.second.size() == 0) { column.first = y; - column.second.resize(1, double(0)); + column.second.resize(1, T(0)); } else if (y < column.first) { column.second.insert(column.second.begin(), - column.first-y, double(0)); + column.first-y, T(0)); column.first = y; } else if (y >= column.first + column.second.size()) { - column.second.resize(y - column.first + 1, double(0)); + column.second.resize(y - column.first + 1, T(0)); } return column.second[y - column.first]; } - SparseMapData &operator+=(const SparseMapData &r); - SparseMapData &operator-=(const SparseMapData &r); - SparseMapData &operator*=(const SparseMapData &r); - SparseMapData &operator/=(const SparseMapData &r); + SparseMapData &operator+=(const SparseMapData &r); + SparseMapData &operator-=(const SparseMapData &r); + SparseMapData &operator*=(const SparseMapData &r); + SparseMapData &operator/=(const SparseMapData &r); - SparseMapData &operator+=(const DenseMapData &r); - SparseMapData &operator-=(const DenseMapData &r); - SparseMapData &operator*=(const DenseMapData &r); - SparseMapData &operator/=(const DenseMapData &r); + SparseMapData &operator+=(const DenseMapData &r); + SparseMapData &operator-=(const DenseMapData &r); + SparseMapData &operator*=(const DenseMapData &r); + SparseMapData &operator/=(const DenseMapData &r); - SparseMapData &operator*=(double r); - SparseMapData &operator/=(double r); + SparseMapData &operator*=(T r); + SparseMapData &operator/=(T r); - SparseMapData *clone(bool copy_data = true) { + SparseMapData *clone(bool copy_data = true) { SparseMapData *m = new SparseMapData(xlen_, ylen_); if (copy_data) { m->data_ = data_; @@ -97,6 +96,18 @@ class SparseMapData { } return m; } + template + SparseMapData *clone_to_type() { + SparseMapData *m = new SparseMapData(xlen_, ylen_); + m->offset_ = offset_; + m->data_.resize(data_.size()); + for (size_t i = 0; i < data_.size(); i++) { + m->data_[i].first = data_[i].first; + m->data_[i].second.resize(data_[i].second.size()); + for (size_t j = 0; j < data_[i].second.size(); j++) + m->data_[i].second[j] = data_[i].second[j]; + } + } DenseMapData *to_dense() const; void compact(); @@ -111,7 +122,7 @@ class SparseMapData { class const_iterator { public: - const_iterator(const SparseMapData &sparse, size_t x_, size_t y_) : + const_iterator(const SparseMapData &sparse, size_t x_, size_t y_) : x(x_), y(y_), sparse_(sparse) {} size_t x, y; @@ -123,13 +134,13 @@ class SparseMapData { return ((x != other.x) || (y != other.y)); } - double operator*() const { return sparse_.at(x, y); } + T operator*() const { return sparse_.at(x, y); } const_iterator operator++(); const_iterator operator++(int) { const_iterator i = *this; ++(*this); return i; } private: - const SparseMapData & sparse_; + const SparseMapData & sparse_; }; const_iterator begin() const { @@ -149,12 +160,11 @@ class SparseMapData { private: uint64_t xlen_, ylen_; - typedef std::pair > data_element; + typedef std::pair > data_element; std::vector data_; uint64_t offset_; }; - class DenseMapData { public: DenseMapData(size_t xlen, size_t ylen) : @@ -212,10 +222,10 @@ class DenseMapData { DenseMapData &operator*=(const DenseMapData &r); DenseMapData &operator/=(const DenseMapData &r); - DenseMapData &operator+=(const SparseMapData &r); - DenseMapData &operator-=(const SparseMapData &r); - DenseMapData &operator*=(const SparseMapData &r); - DenseMapData &operator/=(const SparseMapData &r); + DenseMapData &operator+=(const SparseMapData &r); + DenseMapData &operator-=(const SparseMapData &r); + DenseMapData &operator*=(const SparseMapData &r); + DenseMapData &operator/=(const SparseMapData &r); DenseMapData &operator+=(double r); DenseMapData &operator-=(double r); @@ -288,5 +298,6 @@ class DenseMapData { } }; -CEREAL_CLASS_VERSION(SparseMapData, 1); +CEREAL_CLASS_VERSION(SparseMapData, 1); +CEREAL_CLASS_VERSION(SparseMapData, 1); CEREAL_CLASS_VERSION(DenseMapData, 1); diff --git a/maps/src/maputils.cxx b/maps/src/maputils.cxx index cb075a24..4cd306bd 100644 --- a/maps/src/maputils.cxx +++ b/maps/src/maputils.cxx @@ -19,36 +19,6 @@ using namespace G3Units; -G3SkyMapPtr GetMaskMap(G3SkyMapConstPtr m, bool zero_nans, bool zero_infs) -{ - G3SkyMapPtr mask = m->Clone(false); - - for (size_t i = 0; i < m->size(); i++) { - double v = m->at(i); - if (zero_nans && v != v) - continue; - if (zero_infs && !std::isfinite(v)) - continue; - if (v == 0) - continue; - (*mask)[i] = 1.0; - } - - return mask; -} - -void ApplyMask(G3SkyMapPtr m, G3SkyMapConstPtr mask, bool inverse) -{ - g3_assert(m->IsCompatible(*mask)); - - for (size_t i = 0; i < m->size(); i++) { - if (!(m->at(i))) - continue; - if (!!(mask->at(i)) == inverse) - (*m)[i] = 0; - } -} - void RemoveWeightsT(G3SkyMapPtr T, G3SkyMapWeightsConstPtr W, bool zero_nans) { RemoveWeights(T, NULL, NULL, W, zero_nans); @@ -205,14 +175,10 @@ static double wrap_ra(double ra) } -G3SkyMapPtr GetRaDecMask(G3SkyMapConstPtr m, double ra_left, double ra_right, +G3SkyMapMaskPtr GetRaDecMask(G3SkyMapConstPtr m, double ra_left, double ra_right, double dec_bottom, double dec_top) { - G3SkyMapPtr mask = m->Clone(false); - mask->weighted = false; - mask->units = G3Timestream::None; - mask->pol_type = G3SkyMap::None; - mask->SetPolConv(G3SkyMap::ConvNone); + G3SkyMapMaskPtr mask(new G3SkyMapMask(*m)); ra_left = wrap_ra(ra_left); ra_right = wrap_ra(ra_right); @@ -227,7 +193,7 @@ G3SkyMapPtr GetRaDecMask(G3SkyMapConstPtr m, double ra_left, double ra_right, double dec = radec[1]; if (dec <= dec_bottom || dec >= dec_top) continue; - (*mask)[i] = 1.0; + (*mask)[i] = true; } return mask; @@ -382,7 +348,7 @@ void ReprojMap(G3SkyMapConstPtr in_map, G3SkyMapPtr out_map, int rebin, bool int } // algorithm from https://www.johndcook.com/blog/skewness_kurtosis/ -std::vector GetMapStats(G3SkyMapConstPtr m, G3SkyMapConstPtr mask, +std::vector GetMapMoments(G3SkyMapConstPtr m, G3SkyMapMaskConstPtr mask, int order, bool ignore_zeros, bool ignore_nans, bool ignore_infs) { size_t n = 0; @@ -432,75 +398,8 @@ std::vector GetMapStats(G3SkyMapConstPtr m, G3SkyMapConstPtr mask, } -double -GetMapMedian(G3SkyMapConstPtr m, G3SkyMapConstPtr mask, bool ignore_zeros, bool ignore_nans, - bool ignore_infs) -{ - - std::vector data; - size_t npix = mask ? mask->NpixNonZero() : (ignore_zeros ? m->NpixAllocated() : m->size()); - if (npix == 0) - return 0; - - data.reserve(npix); - - for (size_t i = 0; i < m->size(); i++) { - if (!!mask && !mask->at(i)) - continue; - double v = m->at(i); - if (ignore_zeros && v == 0) - continue; - if (ignore_nans && v != v) - continue; - if (ignore_infs && !std::isfinite(v)) - continue; - data.push_back(v); - } - - npix = data.size(); - size_t n = npix / 2; - std::nth_element(data.begin(), data.begin() + n, data.end()); - - // odd-length array - if (npix % 2) - return data[n]; - - // even-length array - double median = data[n]; - std::nth_element(data.begin(), data.begin() + n - 1, data.end()); - return (median + data[n - 1]) / 2.; -} - - -std::vector -GetMapMinMax(G3SkyMapConstPtr m, G3SkyMapConstPtr mask, bool ignore_zeros, bool ignore_nans, - bool ignore_infs) -{ - double min = 0.0 / 0.0; - double max = 0.0 / 0.0; - - for (size_t i = 0; i < m->size(); i++) { - if (!!mask && !mask->at(i)) - continue; - double v = m->at(i); - if (ignore_zeros && v == 0) - continue; - if (ignore_nans && v != v) - continue; - if (ignore_infs && !std::isfinite(v)) - continue; - if (v > max || max != max) - max = v; - if (v < min || min != min) - min = v; - } - - return {min, max}; -} - - std::vector -GetMapHist(G3SkyMapConstPtr m, const std::vector &bin_edges, G3SkyMapConstPtr mask, +GetMapHist(G3SkyMapConstPtr m, const std::vector &bin_edges, G3SkyMapMaskConstPtr mask, bool ignore_zeros, bool ignore_nans, bool ignore_infs) { @@ -612,37 +511,26 @@ pyconvolve_map(FlatSkyMapConstPtr map, bp::object val) } -void -MakePointSourceMask(G3SkyMapPtr map, const std::vector & ra, +G3SkyMapMaskPtr +MakePointSourceMask(G3SkyMapConstPtr map, const std::vector & ra, const std::vector & dec, const std::vector & radius) { + G3SkyMapMaskPtr mask(new G3SkyMapMask(*map)); g3_assert(ra.size() == dec.size()); g3_assert(ra.size() == radius.size()); for (size_t i = 0; i < ra.size(); i++) { auto pixels = map->QueryDisc(ra[i], dec[i], radius[i]); for (auto p: pixels) - (*map)[p] = 1; + (*mask)[p] = true; } - map->weighted = false; - map->units = G3Timestream::None; - map->pol_type = G3SkyMap::None; - map->SetPolConv(G3SkyMap::ConvNone); + return mask; } namespace bp = boost::python; void maputils_pybindings(void){ - bp::def("get_mask_map", GetMaskMap, - (bp::arg("map_in"), bp::arg("zero_nans")=false, bp::arg("zero_infs")=false), - "Returns a map that is 1 where the input map is nonzero, and optionally " - "non-nan and/or non-infinite"); - - bp::def("apply_mask", ApplyMask, - (bp::arg("map"), bp::arg("mask"), bp::arg("inverse")=false), - "Apply a boolean mask to the input map in place, optionally inverting the mask values"); - bp::def("remove_weights_t", RemoveWeightsT, (bp::arg("T"), bp::arg("W"), bp::arg("zero_nans")=false), "Remove weights from unpolarized maps. If zero_nans is true, empty pixels " @@ -693,8 +581,8 @@ void maputils_pybindings(void){ "polarization conventions. If output attributes are not set, they will be " "copied from the input map."); - bp::def("get_map_stats", GetMapStats, - (bp::arg("map"), bp::arg("mask")=G3SkyMapConstPtr(), bp::arg("order")=2, + bp::def("get_map_moments", GetMapMoments, + (bp::arg("map"), bp::arg("mask")=G3SkyMapMaskConstPtr(), bp::arg("order")=2, bp::arg("ignore_zeros")=false, bp::arg("ignore_nans")=false, bp::arg("ignore_infs")=false), "Computes moment statistics of the input map, optionally ignoring " "zero, nan and/or inf values in the map. If order = 1, only the mean is " @@ -702,22 +590,8 @@ void maputils_pybindings(void){ "are also included, respectively. If a mask is supplied, then only " "the non-zero pixels in the mask are included."); - bp::def("get_map_median", GetMapMedian, - (bp::arg("map"), bp::arg("mask")=G3SkyMapConstPtr(), - bp::arg("ignore_zeros")=false, bp::arg("ignore_nans")=false, bp::arg("ignore_infs")=false), - "Computes the median of the input map, optionally ignoring zero, nan and/or inf " - "values in the map. Requires making a copy of the data. If a mask is " - "supplied, then only the non-zero pixels in the mask are included."); - - bp::def("get_map_minmax", GetMapMinMax, - (bp::arg("map"), bp::arg("mask")=G3SkyMapConstPtr(), - bp::arg("ignore_zeros")=false, bp::arg("ignore_nans")=false, bp::arg("ignore_infs")=false), - "Computes the min and max of the input map, optionally ignoring " - "zero, nan and/or inf values in the map. If a mask is supplied, then " - "only the non-zero pixels in the mask are included."); - bp::def("get_map_hist", GetMapHist, - (bp::arg("map"), bp::arg("bin_edges"), bp::arg("mask")=G3SkyMapConstPtr(), + (bp::arg("map"), bp::arg("bin_edges"), bp::arg("mask")=G3SkyMapMaskConstPtr(), bp::arg("ignore_zeros")=false, bp::arg("ignore_nans")=false, bp::arg("ignore_infs")=false), "Computes the histogram of the input map into bins defined by the array of " "bin edges, optionally ignoring zero, nan and/or inf values in the map. " @@ -729,6 +603,6 @@ void maputils_pybindings(void){ bp::def("make_point_source_mask", MakePointSourceMask, (bp::arg("map"), bp::arg("ra"), bp::arg("dec"), bp::arg("radius")), - "Construct a map with pixels within the given radius around each " - "point source position set to 1."); + "Construct a mask from the input stub map with pixels within the given " + "radius around each point source position set to 1."); } diff --git a/maps/tests/flatsky_maps.py b/maps/tests/flatsky_maps.py index 552c4c4f..33d5358a 100755 --- a/maps/tests/flatsky_maps.py +++ b/maps/tests/flatsky_maps.py @@ -121,6 +121,5 @@ masked |= set(pix2) assert not (set(pix1) ^ set(pix2)) -mask = x.clone(False) -maps.make_point_source_mask(mask, numpy.asarray(avec), numpy.asarray(dvec), numpy.asarray(rvec)) -assert not (set(numpy.where(numpy.asarray(mask).ravel())[0]) ^ masked) +mask = maps.make_point_source_mask(x, numpy.asarray(avec), numpy.asarray(dvec), numpy.asarray(rvec)) +assert not (set(mask.nonzero()) ^ masked) diff --git a/maps/tests/flatskymap_operators.py b/maps/tests/flatskymap_operators.py index cdaec025..1c23f16a 100755 --- a/maps/tests/flatskymap_operators.py +++ b/maps/tests/flatskymap_operators.py @@ -2,7 +2,7 @@ import numpy as np from spt3g import core -from spt3g.maps import FlatSkyMap, MapProjection, get_ra_dec_map, get_map_stats, get_map_median +from spt3g.maps import FlatSkyMap, MapProjection, get_ra_dec_map, get_map_moments from scipy.stats import skew, kurtosis # Sparse extension operators @@ -81,7 +81,8 @@ m *= 2 # Get numbers bigger assert((m == n).all()) assert((m > 0).any()) -assert((m > 0).npix_allocated == 1) +assert((m > 0).sum() == 1) +assert((m > 0).to_map().npix_allocated == 1) m1 = m m2 = m.copy() @@ -252,23 +253,29 @@ # statistics m1 = np.asarray(m).ravel() stats0 = [np.mean(m1), np.var(m1), skew(m1), kurtosis(m1)] - stats1 = get_map_stats(m, order=4) + stats1 = get_map_moments(m, order=4) assert(np.allclose(stats1, stats0)) - med0 = np.median(m1) - med1 = get_map_median(m) + med0 = np.median(np.asarray(m1)) + med1 = np.median(m) assert(np.allclose(med1, med0)) - stats2 = get_map_stats(mpad, order=4, ignore_zeros=True) + stats2 = get_map_moments(mpad, order=4, ignore_zeros=True) assert(np.allclose(stats2, stats0)) - med2 = get_map_median(mpad, ignore_zeros=True) + med2 = mpad.median(where=mpad.to_mask()) assert(np.allclose(med2, med0)) - np.asarray(mpad)[np.asarray(mpad) == 0] = np.nan - stats3 = get_map_stats(mpad, order=4, ignore_nans=True) + mpad[mpad == 0] = np.nan + stats3 = get_map_moments(mpad, order=4, ignore_nans=True) assert(np.allclose(stats3, stats0)) - med3 = get_map_median(mpad, ignore_nans=True) + med3 = mpad.median(where=~mpad.isnan()) assert(np.allclose(med3, med0)) + mask = mpad.to_mask(zero_nans=True) + stats4 = get_map_moments(mpad, mask=mask, order=4) + assert(np.allclose(stats4, stats0)) + med4 = mpad.median(where=mask) + assert(np.allclose(med4, med0)) + # convolution from scipy.signal import convolve2d from spt3g.maps import convolve_map diff --git a/maps/tests/healpix_maps.py b/maps/tests/healpix_maps.py index 3d228785..50f8f48c 100755 --- a/maps/tests/healpix_maps.py +++ b/maps/tests/healpix_maps.py @@ -97,9 +97,8 @@ masked |= set(pix2) assert not (set(pix1) ^ set(pix2)) -mask = x.clone(False) -maps.make_point_source_mask(mask, np.asarray(avec), np.asarray(dvec), np.asarray(rvec)) -assert not (set(np.where(np.asarray(mask).ravel())[0]) ^ masked) +mask = maps.make_point_source_mask(x, np.asarray(avec), np.asarray(dvec), np.asarray(rvec)) +assert not (set(mask.nonzero()) ^ masked) x.shift_ra = True diff --git a/maps/tests/healpixskymap_operators.py b/maps/tests/healpixskymap_operators.py index f0fa1423..a11b3cb2 100755 --- a/maps/tests/healpixskymap_operators.py +++ b/maps/tests/healpixskymap_operators.py @@ -77,7 +77,8 @@ m *= 2 # Get numbers bigger assert((m == n).all()) assert((m > 0).any()) -assert((m > 0).npix_allocated == 1) +assert((m > 0).sum() == 1) +assert((m > 0).to_map().npix_allocated == 1) m1 = m m2 = m.copy() diff --git a/maps/tests/mask_operators.py b/maps/tests/mask_operators.py new file mode 100755 index 00000000..7ede196e --- /dev/null +++ b/maps/tests/mask_operators.py @@ -0,0 +1,190 @@ +#!/usr/bin/env python + +import numpy as np +from spt3g import core, maps + +maplist = [ + maps.FlatSkyMap(np.random.randn(64, 12 * 64), core.G3Units.arcmin), + maps.HealpixSkyMap(np.random.randn(12 * 64 * 64)), +] + +for x in maplist: + print(x) + m = x.to_mask() + + # attributes + assert(m.size == x.size) + assert(m.shape == x.shape) + + # element-wise assignment + if isinstance(x, maps.FlatSkyMap): + v = m[5, 5] + m[5, 5] = not v + assert(m[5, 5] != v) + m[5, 5] = v + assert(m[5, 5] == v) + else: + v = m[5] + m[5] = not v + assert(m[5] != v) + m[5] = v + assert(m[5] == v) + + # masked assignment and retrieval + print("assignment") + pospixels = x[x > 0] + assert(len(pospixels) == (x > 0).sum()) + assert((np.asarray(pospixels) > 0).all()) + + x2 = x.__class__(x) + x2[x2 > 0] = 0 + assert((x2 > 0).sum() == 0) + nzero = (x2 == 0).sum() + x2[x2 == 0] = np.abs(np.random.randn((x2 == 0).sum())) + 1 + assert((x2 > 0).sum() == nzero) + assert((x2 == 0).sum() == 0) + + # init and clone + mx = maps.G3SkyMapMask(x, True) + assert((mx == m).all()) + mn = maps.G3SkyMapMask(x, False) + assert(mn.sum() == 0) + a = np.asarray(x).ravel() + ma = m.array_clone(a) + assert((m == ma).all()) + mp = m.array_clone(a > 0) + assert(((x > 0) == mp).all()) + + # float comparison + print("float comp") + m1 = x == 0 + m2 = x != 0 + assert(m1.sum() + m2.sum() == x.size) + + m1 = x <= 0 + m2 = x > 0 + assert(m1.sum() + m2.sum() == x.size) + + m1 = x < 0 + m2 = x >= 0 + assert(m1.sum() + m2.sum() == x.size) + + # map comparison + print("map comp") + y = 2 * x.clone() + + m1 = y == x + m2 = y != x + assert(m1.sum() + m2.sum() == x.size) + + m1 = y <= x + m2 = y > x + assert(m1.sum() + m2.sum() == x.size) + + m1 = y < x + m2 = y >= x + assert(m1.sum() + m2.sum() == x.size) + + # logic operators + print("logic") + assert((m1 == ~m2).all()) + assert(not (m1 & m2).any()) + assert((m1 | m2).all()) + assert((m1 ^ m2).all()) + + m3 = m1.clone() + m3 |= m2 + assert(m3.all()) + m3 = m1.clone() + m3 &= m2 + assert(not m3.any()) + m3 = m1.clone() + m3 ^= m2 + assert(m3.all()) + m3 = m1.clone() + m3.invert() + assert((m3 == m2).all()) + + # convert to mask + print("conversion") + x2 = x.clone() + x2 *= m1 + assert((x2.to_mask() == m1).all()) + assert(((x * m1).to_mask() == m1).all()) + + x3 = x.clone() + x3.apply_mask(m1, inverse=True) + assert((x3.to_mask() == m2).all()) + + x4 = m1.to_map() + assert(np.sum(x4) == m1.sum()) + assert((x4.to_mask() == m1).all()) + + mpix = m1.nonzero() + xpix, _ = x4.nonzero_pixels() + xpixm = x4.nonzero() + + assert(len(set(mpix) ^ set(xpix)) == 0) + assert(len(set(mpix) ^ set(xpixm)) == 0) + + m1.apply_mask(m2) + assert(not m1.any()) + + # ufuncs + print("ufuncs") + if isinstance(x, maps.FlatSkyMap): + x.sparse = True + else: + x.ringsparse = True + m = x < 0 + + for attr in ["all", "any", "sum", "mean", "var", "std", "min", "max", "argmin", "argmax"]: + for w in [None, m]: + if w is not None and attr in ["argmin", "argmax"]: + continue + kwargs = {"where": w} if w is not None else {} + if attr in ["std", "var"]: + kwargs["ddof"] = 1 + try: + v1 = getattr(np, attr)(x, **kwargs) + except TypeError as e: + # ignore errors with older numpy versions + if w is not None and "unexpected keyword argument 'where'" in str(e): + continue + raise + v2 = getattr(x, attr)(**kwargs) + if w is not None: + kwargs["where"] = np.asarray(w.to_map(), dtype=bool) + if attr in ["min", "max"]: + kwargs["initial"] = np.inf if attr == "min" else -np.inf + v3 = getattr(np, attr)(np.asarray(x.copy()), **kwargs) + print(attr, w, v1, v2, v3) + assert(v1 == v2) + assert(np.isclose(v2, v3)) + if isinstance(x, maps.FlatSkyMap): + assert(x.sparse == True) + else: + assert(x.ringsparse == True) + + x[m.nonzero()[0]] = np.nan + + for attr in ["nansum", "nanmean", "nanvar", "nanstd", "nanmin", "nanmax", "nanargmin", "nanargmax"]: + v1 = getattr(np, attr)(x) + v2 = getattr(x, attr)() + v3 = getattr(np, attr)(np.asarray(x.copy())) + print(attr, v1, v2, v3) + assert(np.isclose(v1, v2)) + assert(np.isclose(v2, v3)) + + x[m.nonzero()[1]] = np.inf + + for attr in ["isinf", "isnan", "isfinite"]: + for w in [None, m]: + print(attr, w) + kwargs = {"where": w} if w is not None else {} + m1 = getattr(x, attr)(**kwargs) + if w is not None: + kwargs["where"] = np.asarray(w.to_map(), dtype=bool) + kwargs["out"] = np.zeros_like(kwargs["where"]) + m2 = maps.G3SkyMapMask(x, getattr(np, attr)(x, **kwargs).ravel()) + assert (m1 == m2).all() diff --git a/maps/tests/weights_operators.py b/maps/tests/weights_operators.py index 6e4bc3c9..06b8cdf3 100755 --- a/maps/tests/weights_operators.py +++ b/maps/tests/weights_operators.py @@ -3,7 +3,7 @@ import numpy as np from spt3g import core from spt3g.maps import G3SkyMapWeights, FlatSkyMap, MapPolType, MapPolConv, MapProjection -from spt3g.maps import get_mask_map, remove_weights, remove_weights_t, apply_weights, apply_weights_t, flatten_pol +from spt3g.maps import remove_weights, remove_weights_t, apply_weights, apply_weights_t, flatten_pol for pol in [True, False]: # allocation @@ -116,9 +116,10 @@ mt.compact(zero_nans=True) assert(mt.npix_allocated == 1) assert(np.allclose(mt[15], np.atleast_1d(vec * 10)[0])) - mask = get_mask_map(mt) + mask = mt.to_mask() assert(mask[15] == 1) - assert(mask.npix_allocated == 1) + assert(mask.sum() == 1) + assert(mask.to_map().npix_allocated == 1) if pol: mq.compact(zero_nans=True)