Support for boolean/mask maps

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)

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<bool>`` as the data storage backend, which uses 64x less memory than a dense map (``std::vector<double>``) 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
 =================

maps/include/maps/FlatSkyMap.h

@@ -11,7 +11,7 @@
 #include <maps/FlatSkyProjection.h>
 
 class DenseMapData;
-class SparseMapData;
+template <typename T> 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<uint64_t> &indices,
 	    std::vector<double> &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<double> *sparse_;
 	uint64_t xpix_, ypix_;
 	bool flat_pol_;
 

maps/include/maps/G3SkyMap.h

@@ -7,6 +7,8 @@
 #include <vector>
 #include <map>
 
+#include <maps/G3SkyMapMask.h>
+
 #include <boost/python.hpp>
 
 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<size_t> AnglesToPixels(const std::vector<double> & alphas,
 	    const std::vector<double> & 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,13 +416,16 @@ class G3SkyMapWeights : public G3FrameObject {
 
 	// Mask
 	G3SkyMapWeights &operator*=(const G3SkyMap &rhs);
+	G3SkyMapWeights &operator*=(const G3SkyMapMask &rhs);
 
 	G3SkyMapPtr Det() const;
 	G3SkyMapPtr Cond() const;
 	boost::shared_ptr<G3SkyMapWeights> Inv() const;
 
 	boost::shared_ptr<G3SkyMapWeights> Rebin(size_t scale) const;
 
+	void ApplyMask(const G3SkyMapMask &mask, bool inverse=false);
+
 	void Compact(bool zero_nans = false);
 
 	boost::shared_ptr<G3SkyMapWeights> Clone(bool copy_data) const {

maps/include/maps/G3SkyMapMask.h

@@ -0,0 +1,128 @@
+#ifndef _G3_SKYMAPMASK_H
+#define _G3_SKYMAPMASK_H
+
+#include <vector>
+
+#include <boost/python.hpp>
+
+/*
+ * 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<G3SkyMapMask> Clone(bool copy_data = true) const;
+	boost::shared_ptr<G3SkyMapMask> ArrayClone(boost::python::object v,
+	    bool zero_nans = false, bool zero_infs = false) const;
+
+	// Return a (modifiable) pixel value
+	std::vector<bool>::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<uint64_t> 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<const G3SkyMap> Parent() const { return parent_; }
+
+	// Return a 1-or-0 sky-map with the contents of the mask (e.g. for plotting)
+	boost::shared_ptr<G3SkyMap> MakeBinaryMap() const;
+
+	class const_iterator {
+	public:
+		typedef std::pair<uint64_t, bool> 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 <class A> void serialize(A &ar, const unsigned v);
+	friend class cereal::access;
+
+	std::vector<bool> data_;
+	boost::shared_ptr<const G3SkyMap> 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
+

maps/include/maps/HealpixSkyMap.h

@@ -10,7 +10,7 @@
 #include <maps/G3SkyMap.h>
 #include <maps/HealpixSkyMapInfo.h>
 
-class SparseMapData;
+template <typename T> 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<uint64_t> &indices,
 	    std::vector<double> &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<double> *dense_;
-	SparseMapData *ring_sparse_;
+	SparseMapData<double> *ring_sparse_;
 	std::unordered_map<uint64_t, double> *indexed_sparse_;
 
 	SET_LOGGER("HealpixSkyMap");