From 42f69c42b014119bdb67fdd42fa8e8bee399b732 Mon Sep 17 00:00:00 2001 From: Nathan Whitehorn Date: Fri, 30 Apr 2021 12:08:51 -0400 Subject: [PATCH 01/55] First, do not harm: template SparseMapData for either bools or doubles. --- maps/include/maps/FlatSkyMap.h | 4 +- maps/include/maps/HealpixSkyMap.h | 4 +- maps/src/FlatSkyMap.cxx | 12 ++--- maps/src/HealpixSkyMap.cxx | 14 +++--- maps/src/mapdata.cxx | 77 +++++++++++++++++++------------ maps/src/mapdata.h | 69 +++++++++++++++------------ 6 files changed, 104 insertions(+), 76 deletions(-) diff --git a/maps/include/maps/FlatSkyMap.h b/maps/include/maps/FlatSkyMap.h index 96368d77..e9b9a77e 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: @@ -186,7 +186,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/HealpixSkyMap.h b/maps/include/maps/HealpixSkyMap.h index 82c566be..7e9ee23a 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 { @@ -140,7 +140,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/src/FlatSkyMap.cxx b/maps/src/FlatSkyMap.cxx index 853ed1b9..d2413ae7 100644 --- a/maps/src/FlatSkyMap.cxx +++ b/maps/src/FlatSkyMap.cxx @@ -80,7 +80,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() @@ -158,7 +158,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_); @@ -282,7 +282,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; @@ -315,7 +315,7 @@ FlatSkyMap::ConvertToSparse() if (!dense_) return; - sparse_ = new SparseMapData(*dense_); + sparse_ = new SparseMapData(*dense_); delete dense_; dense_ = NULL; } @@ -366,7 +366,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); } @@ -410,7 +410,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 \ diff --git a/maps/src/HealpixSkyMap.cxx b/maps/src/HealpixSkyMap.cxx index 3635eb66..2fc0ebf7 100644 --- a/maps/src/HealpixSkyMap.cxx +++ b/maps/src/HealpixSkyMap.cxx @@ -124,7 +124,7 @@ HealpixSkyMap::HealpixSkyMap(boost::python::object v, bool weighted, shift_ra = (phi_max - phi_min) > (phi_max_shift - phi_min_shift); info_.initialize(nside, nested, shift_ra); - ring_sparse_ = new SparseMapData(info_.nring(), info_.nring()); + ring_sparse_ = new SparseMapData(info_.nring(), info_.nring()); for (size_t i = 0; i < sz; i++) (*this)[((unsigned long *)indexview.buf)[i]]= @@ -201,7 +201,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_); @@ -279,7 +279,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: @@ -360,7 +360,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; @@ -405,7 +405,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_; @@ -529,7 +529,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); @@ -970,7 +970,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); 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); From 77754c6178fae4e7bb0864d72db7995dca0f1c1e Mon Sep 17 00:00:00 2001 From: Nathan Whitehorn Date: Tue, 12 Jul 2022 16:40:44 -0400 Subject: [PATCH 02/55] Skymap mask class definition. --- maps/include/maps/G3SkyMapMask.h | 58 ++++++++++++++++++++++++++++++++ 1 file changed, 58 insertions(+) create mode 100644 maps/include/maps/G3SkyMapMask.h diff --git a/maps/include/maps/G3SkyMapMask.h b/maps/include/maps/G3SkyMapMask.h new file mode 100644 index 00000000..86706837 --- /dev/null +++ b/maps/include/maps/G3SkyMapMask.h @@ -0,0 +1,58 @@ +#ifndef _G3_SKYMAP_H +#define _G3_SKYMAP_H + +#include +#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. + */ + +class G3SkyMapMask : public G3FrameObject { +public: + G3SkyMapMask(const G3SkyMap &parent); + virtual ~G3SkyMapMask() {}; + + // Return a (modifiable) pixel value + bool &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 ~= + + G3SkyMapMask operator ~(); + G3SkyMapMask operator |(const G3SkyMapMask &rhs); + G3SkyMapMask operator &(const G3SkyMapMask &rhs); + G3SkyMapMask operator ^(const G3SkyMapMask &rhs); + + bool IsDense() const { return true; } + + // The map for projection info + G3SkyMapConstPtr Parent() const; + +private: + G3SkyMapMask() {} // Fake out for serialization + template void serialize(A &ar, const unsigned v); + friend class cereal::access; + + std::vector data_; + G3SkyMapConstPtr parent_; + + SET_LOGGER("G3SkyMapMask"); +}; + +G3_POINTERS(G3SkyMapMask); +G3_SERIALIZABLE(G3SkyMapMask, 1); + +#endif + From e09f22cdca3920db5000298216ddce90138a087c Mon Sep 17 00:00:00 2001 From: Nathan Whitehorn Date: Tue, 12 Jul 2022 17:21:22 -0400 Subject: [PATCH 03/55] Compilable skymap mask --- maps/CMakeLists.txt | 1 + maps/include/maps/G3SkyMapMask.h | 10 ++-- maps/src/G3SkyMapMask.cxx | 78 ++++++++++++++++++++++++++++++++ 3 files changed, 84 insertions(+), 5 deletions(-) create mode 100644 maps/src/G3SkyMapMask.cxx diff --git a/maps/CMakeLists.txt b/maps/CMakeLists.txt index d0d452c9..8ddf08e5 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 diff --git a/maps/include/maps/G3SkyMapMask.h b/maps/include/maps/G3SkyMapMask.h index 86706837..00ed666b 100644 --- a/maps/include/maps/G3SkyMapMask.h +++ b/maps/include/maps/G3SkyMapMask.h @@ -1,7 +1,7 @@ -#ifndef _G3_SKYMAP_H -#define _G3_SKYMAP_H +#ifndef _G3_SKYMAPMASK_H +#define _G3_SKYMAPMASK_H -#include +#include #include #include @@ -19,7 +19,7 @@ class G3SkyMapMask : public G3FrameObject { virtual ~G3SkyMapMask() {}; // Return a (modifiable) pixel value - bool &operator [] (size_t i); + std::vector::reference operator [] (size_t i); bool operator [] (size_t i) const { return this->at(i); }; bool at(size_t i) const; @@ -38,7 +38,7 @@ class G3SkyMapMask : public G3FrameObject { bool IsDense() const { return true; } // The map for projection info - G3SkyMapConstPtr Parent() const; + G3SkyMapConstPtr Parent() const { return parent_; } private: G3SkyMapMask() {} // Fake out for serialization diff --git a/maps/src/G3SkyMapMask.cxx b/maps/src/G3SkyMapMask.cxx new file mode 100644 index 00000000..7fe9c16d --- /dev/null +++ b/maps/src/G3SkyMapMask.cxx @@ -0,0 +1,78 @@ +#include +#include +#include +#include + +#include +#include + +G3SkyMapMask::G3SkyMapMask(const G3SkyMap &parent) : G3FrameObject() +{ + parent_ = parent.Clone(false); + data_ = std::vector(parent.size()); +} + +std::vector::reference +G3SkyMapMask::operator [] (size_t i) +{ + return data_[i]; +} + +bool +G3SkyMapMask::at (size_t i) const +{ + return data_[i]; +} + +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, int i) +{ + if (i < 0) + i = m.Parent()->size() + i; + if (size_t(i) >= m.Parent()->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, int i, bool val) +{ + if (i < 0) + i = m.Parent()->size() + i; + if (size_t(i) >= m.Parent()->size()) { + PyErr_SetString(PyExc_IndexError, "Index out of range"); + boost::python::throw_error_already_set(); + } + + m[i] = val; +} + +PYBINDINGS("maps") +{ + using namespace boost::python; + + EXPORT_FRAMEOBJECT(G3SkyMapMask, init(), + "Boolean mask of a sky map. Set pixels to use to true, pixels to " + "ignore to false.") + .def_readonly("parent", &G3SkyMapMask::Parent, "\"Parent\" map which " + "contains no data, but can be used to retrieve the parameters of " + "the map to which this mask corresponds.") + .def("__getitem__", &skymapmask_getitem) + .def("__setitem__", &skymapmask_setitem) + ; +} + From 04f52cc89b9a4424ba55e865765df92c41386caa Mon Sep 17 00:00:00 2001 From: Nathan Whitehorn Date: Tue, 12 Jul 2022 17:30:03 -0400 Subject: [PATCH 04/55] Fix copy/paste issue --- maps/src/G3SkyMapMask.cxx | 1 - 1 file changed, 1 deletion(-) diff --git a/maps/src/G3SkyMapMask.cxx b/maps/src/G3SkyMapMask.cxx index 7fe9c16d..743ec1be 100644 --- a/maps/src/G3SkyMapMask.cxx +++ b/maps/src/G3SkyMapMask.cxx @@ -4,7 +4,6 @@ #include #include -#include G3SkyMapMask::G3SkyMapMask(const G3SkyMap &parent) : G3FrameObject() { From 6a1cdc3b6cc62a42612ed4b511006ca1f528d862 Mon Sep 17 00:00:00 2001 From: Nathan Whitehorn Date: Wed, 13 Jul 2022 10:28:57 -0400 Subject: [PATCH 05/55] Finish implementing G3SkyMapMask class. --- maps/include/maps/G3SkyMapMask.h | 4 ++ maps/src/G3SkyMapMask.cxx | 111 +++++++++++++++++++++++++++++++ 2 files changed, 115 insertions(+) diff --git a/maps/include/maps/G3SkyMapMask.h b/maps/include/maps/G3SkyMapMask.h index 00ed666b..5eb2370d 100644 --- a/maps/include/maps/G3SkyMapMask.h +++ b/maps/include/maps/G3SkyMapMask.h @@ -16,6 +16,7 @@ class G3SkyMapMask : public G3FrameObject { public: G3SkyMapMask(const G3SkyMap &parent); + G3SkyMapMask(const G3SkyMapMask &); virtual ~G3SkyMapMask() {}; // Return a (modifiable) pixel value @@ -35,6 +36,9 @@ class G3SkyMapMask : public G3FrameObject { G3SkyMapMask operator &(const G3SkyMapMask &rhs); G3SkyMapMask operator ^(const G3SkyMapMask &rhs); + // Information + bool IsCompatible(const G3SkyMap &map) { return map.IsCompatible(*Parent()); } + bool IsCompatible(const G3SkyMapMask &mask) { return mask.Parent()->IsCompatible(*Parent()); } bool IsDense() const { return true; } // The map for projection info diff --git a/maps/src/G3SkyMapMask.cxx b/maps/src/G3SkyMapMask.cxx index 743ec1be..208c979b 100644 --- a/maps/src/G3SkyMapMask.cxx +++ b/maps/src/G3SkyMapMask.cxx @@ -11,6 +11,12 @@ G3SkyMapMask::G3SkyMapMask(const G3SkyMap &parent) : G3FrameObject() data_ = std::vector(parent.size()); } +G3SkyMapMask::G3SkyMapMask(const G3SkyMapMask &m) : G3FrameObject() +{ + parent_ = m.Parent()->Clone(false); + data_ = std::vector(m.data_); +} + std::vector::reference G3SkyMapMask::operator [] (size_t i) { @@ -23,6 +29,94 @@ G3SkyMapMask::at (size_t i) const return data_[i]; } +G3SkyMapMask & +G3SkyMapMask::operator |=(const G3SkyMapMask &rhs) +{ + g3_assert(IsCompatible(rhs)); + + for (size_t i = 0; i < data_.size(); i++) + (*this)[i] = rhs[i] || (*this)[i]; + + return *this; +} + +G3SkyMapMask & +G3SkyMapMask::operator &=(const G3SkyMapMask &rhs) +{ + g3_assert(IsCompatible(rhs)); + + for (size_t i = 0; i < data_.size(); i++) + (*this)[i] = rhs[i] && (*this)[i]; + + return *this; +} + +G3SkyMapMask & +G3SkyMapMask::operator ^=(const G3SkyMapMask &rhs) +{ + g3_assert(IsCompatible(rhs)); + + for (size_t i = 0; i < data_.size(); i++) + (*this)[i] = rhs[i] ^ (*this)[i]; + + return *this; +} + +G3SkyMapMask & +G3SkyMapMask::invert() +{ + for (size_t i = 0; i < data_.size(); i++) + (*this)[i] = !(*this)[i]; + + return *this; +} + +G3SkyMapMask +G3SkyMapMask::operator ~() +{ + G3SkyMapMask mask(*Parent()); + for (size_t i = 0; i < data_.size(); i++) + mask[i] = !(*this)[i]; + + return mask; +} + +G3SkyMapMask +G3SkyMapMask::operator |(const G3SkyMapMask &rhs) +{ + g3_assert(IsCompatible(rhs)); + + G3SkyMapMask mask(*Parent()); + for (size_t i = 0; i < data_.size(); i++) + mask[i] = (*this)[i] || rhs[i]; + + return mask; +} + +G3SkyMapMask +G3SkyMapMask::operator &(const G3SkyMapMask &rhs) +{ + g3_assert(IsCompatible(rhs)); + + G3SkyMapMask mask(*Parent()); + for (size_t i = 0; i < data_.size(); i++) + mask[i] = (*this)[i] && rhs[i]; + + return mask; +} + +G3SkyMapMask +G3SkyMapMask::operator ^(const G3SkyMapMask &rhs) +{ + g3_assert(IsCompatible(rhs)); + + G3SkyMapMask mask(*Parent()); + for (size_t i = 0; i < data_.size(); i++) + mask[i] = (*this)[i] ^ rhs[i]; + + return mask; +} + template void G3SkyMapMask::serialize(A &ar, unsigned v) { @@ -60,6 +154,15 @@ skymapmask_setitem(G3SkyMapMask &m, int i, bool val) m[i] = val; } +static G3SkyMapMaskPtr +skymapmask_pyinvert(G3SkyMapMaskPtr m) +{ + // Reference-counting problems returning references to Python, + // so use shared pointers. + m->invert(); + return m; +} + PYBINDINGS("maps") { using namespace boost::python; @@ -72,6 +175,14 @@ PYBINDINGS("maps") "the map to which this mask corresponds.") .def("__getitem__", &skymapmask_getitem) .def("__setitem__", &skymapmask_setitem) + .def("invert", &skymapmask_pyinvert, "Invert all elements in 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) ; } From 5266d16786c6fd78ed59d344161ebc7e6a4d7f57 Mon Sep 17 00:00:00 2001 From: Nathan Whitehorn Date: Wed, 13 Jul 2022 10:33:07 -0400 Subject: [PATCH 06/55] Bind all members to Python. --- maps/src/G3SkyMapMask.cxx | 2 ++ 1 file changed, 2 insertions(+) diff --git a/maps/src/G3SkyMapMask.cxx b/maps/src/G3SkyMapMask.cxx index 208c979b..fcf7b59f 100644 --- a/maps/src/G3SkyMapMask.cxx +++ b/maps/src/G3SkyMapMask.cxx @@ -173,6 +173,8 @@ PYBINDINGS("maps") .def_readonly("parent", &G3SkyMapMask::Parent, "\"Parent\" map which " "contains no data, but can be used to retrieve the parameters of " "the map to which this mask corresponds.") + .def("is_compatible", &G3SkyMapMask::IsCompatible, "Returns true if the two masks can be applied to the same map.") + .def("is_compatible", &G3SkyMapMask::IsCompatible, "Returns true if this mask can be applied to the given map.") .def("__getitem__", &skymapmask_getitem) .def("__setitem__", &skymapmask_setitem) .def("invert", &skymapmask_pyinvert, "Invert all elements in mask") From 90f4d8ee5183c63df6880b831c2bb828894c338d Mon Sep 17 00:00:00 2001 From: Nathan Whitehorn Date: Wed, 13 Jul 2022 10:34:06 -0400 Subject: [PATCH 07/55] Sasha doesn't like the verb "to be" --- maps/src/G3SkyMapMask.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/maps/src/G3SkyMapMask.cxx b/maps/src/G3SkyMapMask.cxx index fcf7b59f..31e1d93a 100644 --- a/maps/src/G3SkyMapMask.cxx +++ b/maps/src/G3SkyMapMask.cxx @@ -173,8 +173,8 @@ PYBINDINGS("maps") .def_readonly("parent", &G3SkyMapMask::Parent, "\"Parent\" map which " "contains no data, but can be used to retrieve the parameters of " "the map to which this mask corresponds.") - .def("is_compatible", &G3SkyMapMask::IsCompatible, "Returns true if the two masks can be applied to the same map.") - .def("is_compatible", &G3SkyMapMask::IsCompatible, "Returns true if this mask can be applied to the given map.") + .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("invert", &skymapmask_pyinvert, "Invert all elements in mask") From 2fc4edc2ce769378bd0501427e3b4c9baaf3fb02 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Wed, 13 Jul 2022 09:39:11 -0500 Subject: [PATCH 08/55] const --- maps/include/maps/G3SkyMapMask.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/maps/include/maps/G3SkyMapMask.h b/maps/include/maps/G3SkyMapMask.h index 5eb2370d..106a8dc0 100644 --- a/maps/include/maps/G3SkyMapMask.h +++ b/maps/include/maps/G3SkyMapMask.h @@ -37,8 +37,8 @@ class G3SkyMapMask : public G3FrameObject { G3SkyMapMask operator ^(const G3SkyMapMask &rhs); // Information - bool IsCompatible(const G3SkyMap &map) { return map.IsCompatible(*Parent()); } - bool IsCompatible(const G3SkyMapMask &mask) { return mask.Parent()->IsCompatible(*Parent()); } + bool IsCompatible(const G3SkyMap &map) const { return map.IsCompatible(*Parent()); } + bool IsCompatible(const G3SkyMapMask &mask) const { return mask.Parent()->IsCompatible(*Parent()); } bool IsDense() const { return true; } // The map for projection info From 95dafe654596ab682e1d6ced5d01fbf0c51ddcdc Mon Sep 17 00:00:00 2001 From: Nathan Whitehorn Date: Wed, 13 Jul 2022 10:47:43 -0400 Subject: [PATCH 09/55] Fix Sasha's const commit. --- maps/src/G3SkyMapMask.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/maps/src/G3SkyMapMask.cxx b/maps/src/G3SkyMapMask.cxx index 31e1d93a..b005b992 100644 --- a/maps/src/G3SkyMapMask.cxx +++ b/maps/src/G3SkyMapMask.cxx @@ -173,8 +173,8 @@ PYBINDINGS("maps") .def_readonly("parent", &G3SkyMapMask::Parent, "\"Parent\" map which " "contains no data, but can be used to retrieve the parameters of " "the map to which this mask corresponds.") - .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("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("invert", &skymapmask_pyinvert, "Invert all elements in mask") From 3f9f044a79a4e1bb14e1783622efa16be7ab21da Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Wed, 13 Jul 2022 10:00:54 -0500 Subject: [PATCH 10/55] map * mask operators --- maps/include/maps/G3SkyMap.h | 4 +++ maps/src/G3SkyMap.cxx | 65 ++++++++++++++++++++++++++++++------ 2 files changed, 59 insertions(+), 10 deletions(-) diff --git a/maps/include/maps/G3SkyMap.h b/maps/include/maps/G3SkyMap.h index 6b8935b2..267d8433 100644 --- a/maps/include/maps/G3SkyMap.h +++ b/maps/include/maps/G3SkyMap.h @@ -22,6 +22,8 @@ enum MapCoordReference { * interface for the map maker. */ +class G3SkyMapMask; + class G3SkyMap { public: // Following numerical values are important @@ -95,6 +97,7 @@ class G3SkyMap { // * virtual G3SkyMap &operator*=(const G3SkyMap &rhs); + virtual G3SkyMap &operator*=(const G3SkyMapMask &rhs); virtual G3SkyMap &operator*=(double rhs); // / @@ -360,6 +363,7 @@ class G3SkyMapWeights : public G3FrameObject { // Mask G3SkyMapWeights &operator*=(const G3SkyMap &rhs); + G3SkyMapWeights &operator*=(const G3SkyMapMask &rhs); G3SkyMapPtr Det() const; G3SkyMapPtr Cond() const; diff --git a/maps/src/G3SkyMap.cxx b/maps/src/G3SkyMap.cxx index 681ee94a..1dc0894a 100644 --- a/maps/src/G3SkyMap.cxx +++ b/maps/src/G3SkyMap.cxx @@ -2,6 +2,7 @@ #include #include #include +#include #include #include #include @@ -275,6 +276,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++) @@ -510,6 +522,18 @@ skymap_pynoninplace(multd, *=, double) skymap_pynoninplace(div, /=, G3SkyMap &) skymap_pynoninplace(divd, /=, double) +static G3SkyMapPtr +pyskymap_multm(const G3SkyMap &a, const G3SkyMapMask &b) +{ + G3SkyMapPtr rv = a.Clone(false); + for (size_t i = 0; i < a.size(); i++) { + if (b.at(i) && a.at(i) != 0) + (*rv)[i] = a.at(i); + } + + return rv; +} + static G3SkyMapPtr pyskymap_rsubd(const G3SkyMap &a, const double b) { @@ -588,17 +612,15 @@ pyskymap_pow(const G3SkyMap &a, const G3SkyMap &b) } #define skymap_comp(name, oper) \ -static G3SkyMapPtr \ +static G3SkyMapMaskPtr \ pyskymap_##name(const G3SkyMap &a, const G3SkyMap &b) \ { \ g3_assert(a.IsCompatible(b)); \ g3_assert(a.units == b.units); \ - G3SkyMapPtr rv = a.Clone(false); \ - rv->units = G3Timestream::None; \ - rv->weighted = false; \ + G3SkyMapMaskPtr rv(new G3SkyMapMask(a)); \ for (size_t i = 0; i < a.size(); i++) { \ if (a.at(i) oper b.at(i)) \ - (*rv)[i] = 1; \ + (*rv)[i] = true; \ } \ return rv; \ } @@ -611,15 +633,13 @@ skymap_comp(ge, >=) skymap_comp(gt, >) #define skymap_compd(name, oper) \ -static G3SkyMapPtr \ +static G3SkyMapMaskPtr \ pyskymap_##name(const G3SkyMap &a, const double b) \ { \ - G3SkyMapPtr rv = a.Clone(false); \ - rv->units = G3Timestream::None; \ - rv->weighted = false; \ + G3SkyMapMaskPtr rv(new G3SkyMapMask(a)); \ for (size_t i = 0; i < a.size(); i++) { \ if (a.at(i) oper b) \ - (*rv)[i] = 1; \ + (*rv)[i] = true; \ } \ return rv; \ } @@ -692,6 +712,7 @@ G3SkyMapWeights &G3SkyMapWeights::operator op(rhs_type rhs) \ } skymapweights_inplace(*=, const G3SkyMap &); +skymapweights_inplace(*=, const G3SkyMapMask &); skymapweights_inplace(*=, double); skymapweights_inplace(/=, double); @@ -709,6 +730,25 @@ 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()); + 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; +} + MuellerMatrix MuellerMatrix::Inv() const { MuellerMatrix m; @@ -1040,7 +1080,9 @@ PYBINDINGS("maps") { .def("__sub__", &pyskymap_subd) .def("__rsub__", &pyskymap_rsubd) .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) @@ -1118,7 +1160,10 @@ PYBINDINGS("maps") { .def("__add__", &pyskymapweights_add) .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) From 2e546ca0c0e2972c68f359d08846453f100bf35b Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Wed, 13 Jul 2022 10:23:37 -0500 Subject: [PATCH 11/55] faster mask multiplication in derived classes --- maps/include/maps/FlatSkyMap.h | 1 + maps/include/maps/HealpixSkyMap.h | 1 + maps/src/FlatSkyMap.cxx | 14 ++++++++++++++ maps/src/HealpixSkyMap.cxx | 14 ++++++++++++++ 4 files changed, 30 insertions(+) diff --git a/maps/include/maps/FlatSkyMap.h b/maps/include/maps/FlatSkyMap.h index e9b9a77e..9a5f537b 100644 --- a/maps/include/maps/FlatSkyMap.h +++ b/maps/include/maps/FlatSkyMap.h @@ -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; // / diff --git a/maps/include/maps/HealpixSkyMap.h b/maps/include/maps/HealpixSkyMap.h index 7e9ee23a..e07a88af 100644 --- a/maps/include/maps/HealpixSkyMap.h +++ b/maps/include/maps/HealpixSkyMap.h @@ -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; // / diff --git a/maps/src/FlatSkyMap.cxx b/maps/src/FlatSkyMap.cxx index d2413ae7..58ee8710 100644 --- a/maps/src/FlatSkyMap.cxx +++ b/maps/src/FlatSkyMap.cxx @@ -8,6 +8,7 @@ #include #include +#include #include "mapdata.h" @@ -470,6 +471,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) { diff --git a/maps/src/HealpixSkyMap.cxx b/maps/src/HealpixSkyMap.cxx index 2fc0ebf7..bb756ef3 100644 --- a/maps/src/HealpixSkyMap.cxx +++ b/maps/src/HealpixSkyMap.cxx @@ -9,6 +9,7 @@ #include #include #include +#include #include "mapdata.h" @@ -613,6 +614,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) From 82f425c17cbbe1c5f92f40dbbf029e8757e8f0e0 Mon Sep 17 00:00:00 2001 From: Nathan Whitehorn Date: Wed, 13 Jul 2022 11:18:03 -0400 Subject: [PATCH 12/55] Add 2D indexing when masking flat sky maps. --- maps/src/G3SkyMapMask.cxx | 89 +++++++++++++++++++++++++++++++++++---- 1 file changed, 81 insertions(+), 8 deletions(-) diff --git a/maps/src/G3SkyMapMask.cxx b/maps/src/G3SkyMapMask.cxx index b005b992..4d65ad5b 100644 --- a/maps/src/G3SkyMapMask.cxx +++ b/maps/src/G3SkyMapMask.cxx @@ -4,6 +4,7 @@ #include #include +#include G3SkyMapMask::G3SkyMapMask(const G3SkyMap &parent) : G3FrameObject() { @@ -129,11 +130,47 @@ void G3SkyMapMask::serialize(A &ar, unsigned v) G3_SERIALIZABLE_CODE(G3SkyMapMask); static bool -skymapmask_getitem(const G3SkyMapMask &m, int i) +skymapmask_getitem(const G3SkyMapMask &m, boost::python::object index) { - if (i < 0) - i = m.Parent()->size() + i; - if (size_t(i) >= m.Parent()->size()) { + using namespace boost::python; + + int i = 0; + if (extract(index).check()) { + i = extract(index)(); + + if (i < 0) + i = m.Parent()->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.Parent()->size()) { PyErr_SetString(PyExc_IndexError, "Index out of range"); boost::python::throw_error_already_set(); } @@ -142,11 +179,47 @@ skymapmask_getitem(const G3SkyMapMask &m, int i) } static void -skymapmask_setitem(G3SkyMapMask &m, int i, bool val) +skymapmask_setitem(G3SkyMapMask &m, boost::python::object index, bool val) { - if (i < 0) - i = m.Parent()->size() + i; - if (size_t(i) >= m.Parent()->size()) { + using namespace boost::python; + + int i = 0; + if (extract(index).check()) { + i = extract(index)(); + + if (i < 0) + i = m.Parent()->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.Parent()->size()) { PyErr_SetString(PyExc_IndexError, "Index out of range"); boost::python::throw_error_already_set(); } From aec066f42a982ed4ae6de66135f4766932dc3550 Mon Sep 17 00:00:00 2001 From: Nathan Whitehorn Date: Wed, 13 Jul 2022 12:06:07 -0400 Subject: [PATCH 13/55] Allow initialization of a mask from a map. --- maps/include/maps/G3SkyMapMask.h | 7 ++++++- maps/src/G3SkyMapMask.cxx | 13 ++++++++++--- 2 files changed, 16 insertions(+), 4 deletions(-) diff --git a/maps/include/maps/G3SkyMapMask.h b/maps/include/maps/G3SkyMapMask.h index 106a8dc0..d0e83590 100644 --- a/maps/include/maps/G3SkyMapMask.h +++ b/maps/include/maps/G3SkyMapMask.h @@ -11,11 +11,16 @@ * 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 G3SkyMapMask : public G3FrameObject { public: - G3SkyMapMask(const G3SkyMap &parent); + G3SkyMapMask(const G3SkyMap &parent, bool use_data = false); G3SkyMapMask(const G3SkyMapMask &); virtual ~G3SkyMapMask() {}; diff --git a/maps/src/G3SkyMapMask.cxx b/maps/src/G3SkyMapMask.cxx index 4d65ad5b..c2a32c72 100644 --- a/maps/src/G3SkyMapMask.cxx +++ b/maps/src/G3SkyMapMask.cxx @@ -6,10 +6,15 @@ #include #include -G3SkyMapMask::G3SkyMapMask(const G3SkyMap &parent) : G3FrameObject() +G3SkyMapMask::G3SkyMapMask(const G3SkyMap &parent, bool use_data) + : G3FrameObject() { parent_ = parent.Clone(false); data_ = std::vector(parent.size()); + if (use_data) { + for (size_t i = 0; i < parent.size(); i++) + data_[i] = (parent.at(i) != 0); + } } G3SkyMapMask::G3SkyMapMask(const G3SkyMapMask &m) : G3FrameObject() @@ -240,9 +245,11 @@ PYBINDINGS("maps") { using namespace boost::python; - EXPORT_FRAMEOBJECT(G3SkyMapMask, init(), + EXPORT_FRAMEOBJECT_NOINITNAMESPACE(G3SkyMapMask, (init((arg("parent"), arg("use_data")=false))), "Boolean mask of a sky map. Set pixels to use to true, pixels to " - "ignore to false.") + "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.") .def_readonly("parent", &G3SkyMapMask::Parent, "\"Parent\" map which " "contains no data, but can be used to retrieve the parameters of " "the map to which this mask corresponds.") From 5c9518c6dea3f6445061d1e73b7004f74f09aa83 Mon Sep 17 00:00:00 2001 From: Nathan Whitehorn Date: Wed, 13 Jul 2022 12:14:48 -0400 Subject: [PATCH 14/55] Add conversation to maps. --- maps/include/maps/G3SkyMapMask.h | 3 +++ maps/src/G3SkyMapMask.cxx | 14 ++++++++++++++ 2 files changed, 17 insertions(+) diff --git a/maps/include/maps/G3SkyMapMask.h b/maps/include/maps/G3SkyMapMask.h index d0e83590..1f7d052a 100644 --- a/maps/include/maps/G3SkyMapMask.h +++ b/maps/include/maps/G3SkyMapMask.h @@ -49,6 +49,9 @@ class G3SkyMapMask : public G3FrameObject { // The map for projection info G3SkyMapConstPtr Parent() const { return parent_; } + // Return a 1-or-0 sky-map with the contents of the mask (e.g. for plotting) + G3SkyMapPtr MakeBinaryMap(); + private: G3SkyMapMask() {} // Fake out for serialization template void serialize(A &ar, const unsigned v); diff --git a/maps/src/G3SkyMapMask.cxx b/maps/src/G3SkyMapMask.cxx index c2a32c72..cd7d2268 100644 --- a/maps/src/G3SkyMapMask.cxx +++ b/maps/src/G3SkyMapMask.cxx @@ -123,6 +123,19 @@ G3SkyMapMask::operator ^(const G3SkyMapMask &rhs) return mask; } +G3SkyMapPtr +G3SkyMapMask::MakeBinaryMap() +{ + G3SkyMapPtr out = parent_->Clone(); + + for (size_t i = 0; i < data_.size(); i++) { + if (at(i)) + (*out)[i] = 1.0; + } + + return out; +} + template void G3SkyMapMask::serialize(A &ar, unsigned v) { @@ -265,6 +278,7 @@ PYBINDINGS("maps") .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.") ; } From 45d1919d9198be74f6fa512caf924c22cdc4ebc9 Mon Sep 17 00:00:00 2001 From: Nathan Whitehorn Date: Wed, 13 Jul 2022 12:19:42 -0400 Subject: [PATCH 15/55] More const. --- maps/include/maps/G3SkyMapMask.h | 10 +++++----- maps/src/G3SkyMapMask.cxx | 10 +++++----- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/maps/include/maps/G3SkyMapMask.h b/maps/include/maps/G3SkyMapMask.h index 1f7d052a..d1154b35 100644 --- a/maps/include/maps/G3SkyMapMask.h +++ b/maps/include/maps/G3SkyMapMask.h @@ -36,10 +36,10 @@ class G3SkyMapMask : public G3FrameObject { G3SkyMapMask &operator ^=(const G3SkyMapMask &rhs); G3SkyMapMask &invert(); // Basically ~= - G3SkyMapMask operator ~(); - G3SkyMapMask operator |(const G3SkyMapMask &rhs); - G3SkyMapMask operator &(const G3SkyMapMask &rhs); - G3SkyMapMask operator ^(const G3SkyMapMask &rhs); + G3SkyMapMask operator ~() const; + G3SkyMapMask operator |(const G3SkyMapMask &rhs) const; + G3SkyMapMask operator &(const G3SkyMapMask &rhs) const; + G3SkyMapMask operator ^(const G3SkyMapMask &rhs) const; // Information bool IsCompatible(const G3SkyMap &map) const { return map.IsCompatible(*Parent()); } @@ -50,7 +50,7 @@ class G3SkyMapMask : public G3FrameObject { G3SkyMapConstPtr Parent() const { return parent_; } // Return a 1-or-0 sky-map with the contents of the mask (e.g. for plotting) - G3SkyMapPtr MakeBinaryMap(); + G3SkyMapPtr MakeBinaryMap() const; private: G3SkyMapMask() {} // Fake out for serialization diff --git a/maps/src/G3SkyMapMask.cxx b/maps/src/G3SkyMapMask.cxx index cd7d2268..70dfa6c1 100644 --- a/maps/src/G3SkyMapMask.cxx +++ b/maps/src/G3SkyMapMask.cxx @@ -78,7 +78,7 @@ G3SkyMapMask::invert() } G3SkyMapMask -G3SkyMapMask::operator ~() +G3SkyMapMask::operator ~() const { G3SkyMapMask mask(*Parent()); for (size_t i = 0; i < data_.size(); i++) @@ -88,7 +88,7 @@ G3SkyMapMask::operator ~() } G3SkyMapMask -G3SkyMapMask::operator |(const G3SkyMapMask &rhs) +G3SkyMapMask::operator |(const G3SkyMapMask &rhs) const { g3_assert(IsCompatible(rhs)); @@ -100,7 +100,7 @@ G3SkyMapMask::operator |(const G3SkyMapMask &rhs) } G3SkyMapMask -G3SkyMapMask::operator &(const G3SkyMapMask &rhs) +G3SkyMapMask::operator &(const G3SkyMapMask &rhs) const { g3_assert(IsCompatible(rhs)); @@ -112,7 +112,7 @@ G3SkyMapMask::operator &(const G3SkyMapMask &rhs) } G3SkyMapMask -G3SkyMapMask::operator ^(const G3SkyMapMask &rhs) +G3SkyMapMask::operator ^(const G3SkyMapMask &rhs) const { g3_assert(IsCompatible(rhs)); @@ -124,7 +124,7 @@ G3SkyMapMask::operator ^(const G3SkyMapMask &rhs) } G3SkyMapPtr -G3SkyMapMask::MakeBinaryMap() +G3SkyMapMask::MakeBinaryMap() const { G3SkyMapPtr out = parent_->Clone(); From 333e4e8b1bd615cecd96f676ce729af1e0f84a8f Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Wed, 13 Jul 2022 11:39:28 -0500 Subject: [PATCH 16/55] all, any, sum, nonzero_pixels --- maps/include/maps/G3SkyMapMask.h | 4 ++ maps/src/G3SkyMapMask.cxx | 59 +++++++++++++++++++++++---- maps/tests/flatskymap_operators.py | 3 +- maps/tests/healpixskymap_operators.py | 3 +- 4 files changed, 60 insertions(+), 9 deletions(-) diff --git a/maps/include/maps/G3SkyMapMask.h b/maps/include/maps/G3SkyMapMask.h index d1154b35..77dcd1da 100644 --- a/maps/include/maps/G3SkyMapMask.h +++ b/maps/include/maps/G3SkyMapMask.h @@ -35,6 +35,10 @@ class G3SkyMapMask : public G3FrameObject { G3SkyMapMask &operator &=(const G3SkyMapMask &rhs); G3SkyMapMask &operator ^=(const G3SkyMapMask &rhs); G3SkyMapMask &invert(); // Basically ~= + bool all() const; + bool any() const; + size_t sum() const; + std::vector NonZeroPixels() const; G3SkyMapMask operator ~() const; G3SkyMapMask operator |(const G3SkyMapMask &rhs) const; diff --git a/maps/src/G3SkyMapMask.cxx b/maps/src/G3SkyMapMask.cxx index 70dfa6c1..e9b9df84 100644 --- a/maps/src/G3SkyMapMask.cxx +++ b/maps/src/G3SkyMapMask.cxx @@ -41,7 +41,7 @@ G3SkyMapMask::operator |=(const G3SkyMapMask &rhs) g3_assert(IsCompatible(rhs)); for (size_t i = 0; i < data_.size(); i++) - (*this)[i] = rhs[i] || (*this)[i]; + (*this)[i] = rhs.at(i) || (*this)[i]; return *this; } @@ -52,7 +52,7 @@ G3SkyMapMask::operator &=(const G3SkyMapMask &rhs) g3_assert(IsCompatible(rhs)); for (size_t i = 0; i < data_.size(); i++) - (*this)[i] = rhs[i] && (*this)[i]; + (*this)[i] = rhs.at(i) && (*this)[i]; return *this; } @@ -63,7 +63,7 @@ G3SkyMapMask::operator ^=(const G3SkyMapMask &rhs) g3_assert(IsCompatible(rhs)); for (size_t i = 0; i < data_.size(); i++) - (*this)[i] = rhs[i] ^ (*this)[i]; + (*this)[i] = rhs.at(i) ^ (*this)[i]; return *this; } @@ -77,12 +77,52 @@ G3SkyMapMask::invert() return *this; } +bool +G3SkyMapMask::all() const +{ + for (size_t i = 0; i < data_.size(); i++) + if (at(i) == 0) + return false; + return true; +} + +bool +G3SkyMapMask::any() const +{ + for (size_t i = 0; i < data_.size(); i++) + if (at(i) != 0) + return true; + return false; +} + +size_t +G3SkyMapMask::sum() const +{ + size_t s = 0; + for (size_t i = 0; i < data_.size(); i++) + s += at(i); + return s; +} + +std::vector +G3SkyMapMask::NonZeroPixels() const +{ + std::vector indices; + + for (size_t i = 0; i < data_.size(); i++) { + if (at(i)) + indices.push_back(i); + } + + return indices; +} + G3SkyMapMask G3SkyMapMask::operator ~() const { G3SkyMapMask mask(*Parent()); for (size_t i = 0; i < data_.size(); i++) - mask[i] = !(*this)[i]; + mask[i] = !at(i); return mask; } @@ -94,7 +134,7 @@ G3SkyMapMask::operator |(const G3SkyMapMask &rhs) const G3SkyMapMask mask(*Parent()); for (size_t i = 0; i < data_.size(); i++) - mask[i] = (*this)[i] || rhs[i]; + mask[i] = at(i) || rhs.at(i); return mask; } @@ -106,7 +146,7 @@ G3SkyMapMask::operator &(const G3SkyMapMask &rhs) const G3SkyMapMask mask(*Parent()); for (size_t i = 0; i < data_.size(); i++) - mask[i] = (*this)[i] && rhs[i]; + mask[i] = at(i) && rhs.at(i); return mask; } @@ -118,7 +158,7 @@ G3SkyMapMask::operator ^(const G3SkyMapMask &rhs) const G3SkyMapMask mask(*Parent()); for (size_t i = 0; i < data_.size(); i++) - mask[i] = (*this)[i] ^ rhs[i]; + mask[i] = at(i) ^ rhs.at(i); return mask; } @@ -271,6 +311,11 @@ PYBINDINGS("maps") .def("__getitem__", &skymapmask_getitem) .def("__setitem__", &skymapmask_setitem) .def("invert", &skymapmask_pyinvert, "Invert all elements in mask") + .def("all", &G3SkyMapMask::all, "Test whether all elements are non-zero") + .def("any", &G3SkyMapMask::any, "Test whether any elements are non-zero") + .def("sum", &G3SkyMapMask::sum, "Sum of all elements in mask") + .def("nonzero_pixels", &G3SkyMapMask::NonZeroPixels, + "Return a list of indices of non-zero pixels in the mask") .def(bp::self |= bp::self) .def(bp::self &= bp::self) .def(bp::self ^= bp::self) diff --git a/maps/tests/flatskymap_operators.py b/maps/tests/flatskymap_operators.py index cdaec025..d92b4f3d 100755 --- a/maps/tests/flatskymap_operators.py +++ b/maps/tests/flatskymap_operators.py @@ -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() 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() From 25081c9a4c324290f1a219855b64878b6138ab59 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Wed, 13 Jul 2022 12:11:52 -0500 Subject: [PATCH 17/55] in-place mask multiplication --- CMakeLists.txt | 5 +++++ maps/src/G3SkyMap.cxx | 31 ++++++++++++++++++++++++++++++- 2 files changed, 35 insertions(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 53502668..baf1cb97 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -143,6 +143,11 @@ foreach(subdir ${SUBDIRS}) add_subdirectory(${CMAKE_SOURCE_DIR}/${pname}) endforeach(subdir ${SUBDIRS}) +if ("${CMAKE_BUILD_TYPE}" STREQUAL "Debug") + add_executable(cppexample ${CMAKE_SOURCE_DIR}/examples/cppexample.cxx) + target_link_libraries(cppexample ${SPT3G_LIBRARIES}) +endif("${CMAKE_BUILD_TYPE}" STREQUAL "Debug") + # export configuration files for use in other projects export(TARGETS ${SPT3G_LIBRARIES} NAMESPACE spt3g:: FILE ${CMAKE_BINARY_DIR}/cmake/Spt3gTargets.cmake) configure_file(${CMAKE_SOURCE_DIR}/cmake/Spt3gConfig.cmake.in ${CMAKE_BINARY_DIR}/cmake/Spt3gConfig.cmake @ONLY) diff --git a/maps/src/G3SkyMap.cxx b/maps/src/G3SkyMap.cxx index 1dc0894a..e49e8d60 100644 --- a/maps/src/G3SkyMap.cxx +++ b/maps/src/G3SkyMap.cxx @@ -534,6 +534,13 @@ pyskymap_multm(const G3SkyMap &a, const G3SkyMapMask &b) 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) { @@ -734,7 +741,8 @@ static G3SkyMapWeightsPtr pyskymapweights_multma(const G3SkyMapWeights &a, const G3SkyMapMask &b) { G3SkyMapWeightsPtr rv(new G3SkyMapWeights()); - rv->TT = pyskymap_multm(*a.TT, b); + if (!!a.TT) + rv->TT = pyskymap_multm(*a.TT, b); if (!!a.TQ) rv->TQ = pyskymap_multm(*a.TQ, b); if (!!a.TU) @@ -749,6 +757,25 @@ pyskymapweights_multma(const G3SkyMapWeights &a, const G3SkyMapMask &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; @@ -1079,6 +1106,7 @@ 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) @@ -1159,6 +1187,7 @@ 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) From 7fdfab77305c7fa99b09f12a64e96cc109c76833 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Wed, 13 Jul 2022 12:17:26 -0500 Subject: [PATCH 18/55] revert accidental commit --- CMakeLists.txt | 5 ----- 1 file changed, 5 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index baf1cb97..53502668 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -143,11 +143,6 @@ foreach(subdir ${SUBDIRS}) add_subdirectory(${CMAKE_SOURCE_DIR}/${pname}) endforeach(subdir ${SUBDIRS}) -if ("${CMAKE_BUILD_TYPE}" STREQUAL "Debug") - add_executable(cppexample ${CMAKE_SOURCE_DIR}/examples/cppexample.cxx) - target_link_libraries(cppexample ${SPT3G_LIBRARIES}) -endif("${CMAKE_BUILD_TYPE}" STREQUAL "Debug") - # export configuration files for use in other projects export(TARGETS ${SPT3G_LIBRARIES} NAMESPACE spt3g:: FILE ${CMAKE_BINARY_DIR}/cmake/Spt3gTargets.cmake) configure_file(${CMAKE_SOURCE_DIR}/cmake/Spt3gConfig.cmake.in ${CMAKE_BINARY_DIR}/cmake/Spt3gConfig.cmake @ONLY) From e874b960f3c83296c85cafcab78436e30b3fd863 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Wed, 13 Jul 2022 15:00:23 -0500 Subject: [PATCH 19/55] avoid circular header dependency --- maps/include/maps/G3SkyMap.h | 2 ++ maps/include/maps/G3SkyMapMask.h | 13 +++++++------ maps/src/G3SkyMap.cxx | 11 +++++++++++ maps/src/G3SkyMapMask.cxx | 17 +++++++++++++++-- 4 files changed, 35 insertions(+), 8 deletions(-) diff --git a/maps/include/maps/G3SkyMap.h b/maps/include/maps/G3SkyMap.h index 267d8433..2ba1e03b 100644 --- a/maps/include/maps/G3SkyMap.h +++ b/maps/include/maps/G3SkyMap.h @@ -151,6 +151,8 @@ class G3SkyMap { throw std::runtime_error("Compactification not implemented"); } + virtual boost::shared_ptr MakeMask() const; + protected: MapPolConv pol_conv_; diff --git a/maps/include/maps/G3SkyMapMask.h b/maps/include/maps/G3SkyMapMask.h index 77dcd1da..565d781a 100644 --- a/maps/include/maps/G3SkyMapMask.h +++ b/maps/include/maps/G3SkyMapMask.h @@ -1,7 +1,6 @@ #ifndef _G3_SKYMAPMASK_H #define _G3_SKYMAPMASK_H -#include #include #include @@ -18,6 +17,8 @@ * non-zero. */ +class G3SkyMap; + class G3SkyMapMask : public G3FrameObject { public: G3SkyMapMask(const G3SkyMap &parent, bool use_data = false); @@ -46,15 +47,15 @@ class G3SkyMapMask : public G3FrameObject { G3SkyMapMask operator ^(const G3SkyMapMask &rhs) const; // Information - bool IsCompatible(const G3SkyMap &map) const { return map.IsCompatible(*Parent()); } - bool IsCompatible(const G3SkyMapMask &mask) const { return mask.Parent()->IsCompatible(*Parent()); } + bool IsCompatible(const G3SkyMap &map) const; + bool IsCompatible(const G3SkyMapMask &mask) const; bool IsDense() const { return true; } // The map for projection info - G3SkyMapConstPtr Parent() const { return parent_; } + boost::shared_ptr Parent() const { return parent_; } // Return a 1-or-0 sky-map with the contents of the mask (e.g. for plotting) - G3SkyMapPtr MakeBinaryMap() const; + boost::shared_ptr MakeBinaryMap() const; private: G3SkyMapMask() {} // Fake out for serialization @@ -62,7 +63,7 @@ class G3SkyMapMask : public G3FrameObject { friend class cereal::access; std::vector data_; - G3SkyMapConstPtr parent_; + boost::shared_ptr parent_; SET_LOGGER("G3SkyMapMask"); }; diff --git a/maps/src/G3SkyMap.cxx b/maps/src/G3SkyMap.cxx index e49e8d60..c954cd36 100644 --- a/maps/src/G3SkyMap.cxx +++ b/maps/src/G3SkyMap.cxx @@ -227,6 +227,13 @@ size_t G3SkyMap::size() const return s; } +G3SkyMapMaskPtr +G3SkyMap::MakeMask() const +{ + G3SkyMapMaskPtr m(new G3SkyMapMask(*this, true)); + return m; +} + G3SkyMap &G3SkyMap::operator+=(const G3SkyMap & rhs) { g3_assert(IsCompatible(rhs)); @@ -1091,6 +1098,10 @@ PYBINDINGS("maps") { "that is already sparse will be compactified in place in its " "current representation without additional memory overhead.") + .def("to_mask", &G3SkyMap::MakeMask, + "Create a G3SkyMapMask object from the parent map, with pixels " + "set to true where the map is non-zero.") + .def(bp::self += bp::self) .def(bp::self *= bp::self) .def(bp::self -= bp::self) diff --git a/maps/src/G3SkyMapMask.cxx b/maps/src/G3SkyMapMask.cxx index e9b9df84..09a7d2db 100644 --- a/maps/src/G3SkyMapMask.cxx +++ b/maps/src/G3SkyMapMask.cxx @@ -4,6 +4,7 @@ #include #include +#include #include G3SkyMapMask::G3SkyMapMask(const G3SkyMap &parent, bool use_data) @@ -19,7 +20,7 @@ G3SkyMapMask::G3SkyMapMask(const G3SkyMap &parent, bool use_data) G3SkyMapMask::G3SkyMapMask(const G3SkyMapMask &m) : G3FrameObject() { - parent_ = m.Parent()->Clone(false); + parent_ = m.parent_; data_ = std::vector(m.data_); } @@ -163,10 +164,22 @@ G3SkyMapMask::operator ^(const G3SkyMapMask &rhs) const return mask; } +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(); + G3SkyMapPtr out = Parent()->Clone(); for (size_t i = 0; i < data_.size(); i++) { if (at(i)) From f597cb45041fa97363c398953de20d2289cf370a Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Wed, 13 Jul 2022 16:18:43 -0500 Subject: [PATCH 20/55] comparison operator members of G3SkyMap --- maps/include/maps/G3SkyMap.h | 27 ++++++++++++- maps/src/G3SkyMap.cxx | 73 ++++++++++++++++++------------------ 2 files changed, 62 insertions(+), 38 deletions(-) diff --git a/maps/include/maps/G3SkyMap.h b/maps/include/maps/G3SkyMap.h index 2ba1e03b..d948cac5 100644 --- a/maps/include/maps/G3SkyMap.h +++ b/maps/include/maps/G3SkyMap.h @@ -7,6 +7,8 @@ #include #include +#include + #include enum MapCoordReference { @@ -22,8 +24,6 @@ enum MapCoordReference { * interface for the map maker. */ -class G3SkyMapMask; - class G3SkyMap { public: // Following numerical values are important @@ -104,6 +104,29 @@ class G3SkyMap { 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() const; + virtual bool any() const; + // Pointing information std::vector AnglesToPixels(const std::vector & alphas, const std::vector & deltas) const; diff --git a/maps/src/G3SkyMap.cxx b/maps/src/G3SkyMap.cxx index c954cd36..76f04987 100644 --- a/maps/src/G3SkyMap.cxx +++ b/maps/src/G3SkyMap.cxx @@ -626,15 +626,15 @@ pyskymap_pow(const G3SkyMap &a, const G3SkyMap &b) } #define skymap_comp(name, oper) \ -static G3SkyMapMaskPtr \ -pyskymap_##name(const G3SkyMap &a, const G3SkyMap &b) \ +G3SkyMapMask \ +G3SkyMap::operator oper(const G3SkyMap &rhs) \ { \ - g3_assert(a.IsCompatible(b)); \ - g3_assert(a.units == b.units); \ - G3SkyMapMaskPtr rv(new G3SkyMapMask(a)); \ - for (size_t i = 0; i < a.size(); i++) { \ - if (a.at(i) oper b.at(i)) \ - (*rv)[i] = true; \ + 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; \ } @@ -647,13 +647,13 @@ skymap_comp(ge, >=) skymap_comp(gt, >) #define skymap_compd(name, oper) \ -static G3SkyMapMaskPtr \ -pyskymap_##name(const G3SkyMap &a, const double b) \ +G3SkyMapMask \ +G3SkyMap::operator oper(double rhs) \ { \ - G3SkyMapMaskPtr rv(new G3SkyMapMask(a)); \ - for (size_t i = 0; i < a.size(); i++) { \ - if (a.at(i) oper b) \ - (*rv)[i] = true; \ + G3SkyMapMask rv(*this); \ + for (size_t i = 0; i < size(); i++) { \ + if (at(i) oper rhs) \ + rv[i] = true; \ } \ return rv; \ } @@ -665,20 +665,20 @@ skymap_compd(ned, !=) skymap_compd(ged, >=) skymap_compd(gtd, >) -static bool -pyskymap_all(const G3SkyMap &a) +bool +G3SkyMap::all() const { - for (size_t i = 0; i < a.size(); i++) - if (a.at(i) == 0) + for (size_t i = 0; i < size(); i++) + if (at(i) == 0) return false; return true; } -static bool -pyskymap_any(const G3SkyMap &a) +bool +G3SkyMap::any() const { - for (size_t i = 0; i < a.size(); i++) - if (a.at(i) != 0) + for (size_t i = 0; i < size(); i++) + if (at(i) != 0) return true; return false; } @@ -1134,21 +1134,22 @@ 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("any", &G3SkyMap::any) + .def("all", &G3SkyMap::all) ; EXPORT_FRAMEOBJECT(G3SkyMapWeights, init<>(), From f7cb19c539960f384f02596194012c281fad239c Mon Sep 17 00:00:00 2001 From: Nathan Whitehorn Date: Thu, 14 Jul 2022 10:42:23 -0400 Subject: [PATCH 21/55] Missed key pointer conversions in py bindings. --- maps/src/G3SkyMapMask.cxx | 1 + 1 file changed, 1 insertion(+) diff --git a/maps/src/G3SkyMapMask.cxx b/maps/src/G3SkyMapMask.cxx index 09a7d2db..8abd6126 100644 --- a/maps/src/G3SkyMapMask.cxx +++ b/maps/src/G3SkyMapMask.cxx @@ -338,5 +338,6 @@ PYBINDINGS("maps") .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_ponter_conversions(); } From 6d7fd75b102d48dda91d932f84542de47c30790e Mon Sep 17 00:00:00 2001 From: Nathan Whitehorn Date: Thu, 14 Jul 2022 10:44:07 -0400 Subject: [PATCH 22/55] Fix typo. --- maps/src/G3SkyMapMask.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/maps/src/G3SkyMapMask.cxx b/maps/src/G3SkyMapMask.cxx index 8abd6126..8eccc90f 100644 --- a/maps/src/G3SkyMapMask.cxx +++ b/maps/src/G3SkyMapMask.cxx @@ -338,6 +338,6 @@ PYBINDINGS("maps") .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_ponter_conversions(); + register_pointer_conversions(); } From f62c3d2c1807c671b4b4b463bbff42e3bd3b3cd7 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Thu, 14 Jul 2022 09:20:47 -0500 Subject: [PATCH 23/55] optionally exclude nans and infs in mask constructor --- maps/include/maps/G3SkyMap.h | 3 ++- maps/include/maps/G3SkyMapMask.h | 3 ++- maps/src/G3SkyMap.cxx | 8 +++++--- maps/src/G3SkyMapMask.cxx | 19 +++++++++++++------ 4 files changed, 22 insertions(+), 11 deletions(-) diff --git a/maps/include/maps/G3SkyMap.h b/maps/include/maps/G3SkyMap.h index d948cac5..958cfe35 100644 --- a/maps/include/maps/G3SkyMap.h +++ b/maps/include/maps/G3SkyMap.h @@ -174,7 +174,8 @@ class G3SkyMap { throw std::runtime_error("Compactification not implemented"); } - virtual boost::shared_ptr MakeMask() const; + virtual boost::shared_ptr MakeMask(bool zero_nans = false, + bool zero_infs = false) const; protected: MapPolConv pol_conv_; diff --git a/maps/include/maps/G3SkyMapMask.h b/maps/include/maps/G3SkyMapMask.h index 565d781a..b6564383 100644 --- a/maps/include/maps/G3SkyMapMask.h +++ b/maps/include/maps/G3SkyMapMask.h @@ -21,7 +21,8 @@ class G3SkyMap; class G3SkyMapMask : public G3FrameObject { public: - G3SkyMapMask(const G3SkyMap &parent, bool use_data = false); + G3SkyMapMask(const G3SkyMap &parent, bool use_data = false, + bool zero_nans = false, bool zero_infs = false); G3SkyMapMask(const G3SkyMapMask &); virtual ~G3SkyMapMask() {}; diff --git a/maps/src/G3SkyMap.cxx b/maps/src/G3SkyMap.cxx index 76f04987..c7ae6de4 100644 --- a/maps/src/G3SkyMap.cxx +++ b/maps/src/G3SkyMap.cxx @@ -228,9 +228,9 @@ size_t G3SkyMap::size() const } G3SkyMapMaskPtr -G3SkyMap::MakeMask() const +G3SkyMap::MakeMask(bool zero_nans, bool zero_infs) const { - G3SkyMapMaskPtr m(new G3SkyMapMask(*this, true)); + G3SkyMapMaskPtr m(new G3SkyMapMask(*this, true, zero_nans, zero_infs)); return m; } @@ -1099,8 +1099,10 @@ PYBINDINGS("maps") { "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.") + "set to true where the map is non-zero (and optionally non-nan " + "and/or finite).") .def(bp::self += bp::self) .def(bp::self *= bp::self) diff --git a/maps/src/G3SkyMapMask.cxx b/maps/src/G3SkyMapMask.cxx index 8eccc90f..a798f298 100644 --- a/maps/src/G3SkyMapMask.cxx +++ b/maps/src/G3SkyMapMask.cxx @@ -7,14 +7,20 @@ #include #include -G3SkyMapMask::G3SkyMapMask(const G3SkyMap &parent, bool use_data) - : G3FrameObject() +G3SkyMapMask::G3SkyMapMask(const G3SkyMap &parent, bool use_data, + bool zero_nans, bool zero_infs) : G3FrameObject() { parent_ = parent.Clone(false); data_ = std::vector(parent.size()); if (use_data) { - for (size_t i = 0; i < parent.size(); i++) - data_[i] = (parent.at(i) != 0); + for (size_t i = 0; i < parent.size(); i++) { + double v = parent.at(i); + if (zero_nans && v != v) + continue; + if (zero_infs && !std::isfinite(v)) + continue; + data_[i] = (v != 0); + } } } @@ -311,11 +317,12 @@ PYBINDINGS("maps") { using namespace boost::python; - EXPORT_FRAMEOBJECT_NOINITNAMESPACE(G3SkyMapMask, (init((arg("parent"), arg("use_data")=false))), + EXPORT_FRAMEOBJECT_NOINITNAMESPACE(G3SkyMapMask, (init((arg("parent"), arg("use_data")=false, arg("zero_nans")=false, arg("zero_infs")=false))), "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.") + "initialized to zero. Use zero_nans or zero_infs to exclude nan " + "or inf elements from the mask.") .def_readonly("parent", &G3SkyMapMask::Parent, "\"Parent\" map which " "contains no data, but can be used to retrieve the parameters of " "the map to which this mask corresponds.") From 3bbc22ec13f9d28644f16077d5560f70bcd694fe Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Thu, 14 Jul 2022 10:42:57 -0500 Subject: [PATCH 24/55] use masks in maputils --- maps/include/maps/maputils.h | 22 ++++----- maps/python/map_modules.py | 2 +- maps/src/maputils.cxx | 83 ++++++++------------------------- maps/tests/flatsky_maps.py | 5 +- maps/tests/healpix_maps.py | 5 +- maps/tests/weights_operators.py | 7 +-- 6 files changed, 36 insertions(+), 88 deletions(-) diff --git a/maps/include/maps/maputils.h b/maps/include/maps/maputils.h index ee93ce52..8ce46e85 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,32 @@ 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, +std::vector GetMapStats(G3SkyMapConstPtr m, G3SkyMapMaskConstPtr 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, +double GetMapMedian(G3SkyMapConstPtr m, G3SkyMapMaskConstPtr 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 GetMapMinMax(G3SkyMapConstPtr m, G3SkyMapMaskConstPtr mask=NULL, 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 6c896ce6..69e9b861 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/src/maputils.cxx b/maps/src/maputils.cxx index cb075a24..620c2ad9 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 GetMapStats(G3SkyMapConstPtr m, G3SkyMapMaskConstPtr mask, int order, bool ignore_zeros, bool ignore_nans, bool ignore_infs) { size_t n = 0; @@ -433,12 +399,12 @@ std::vector GetMapStats(G3SkyMapConstPtr m, G3SkyMapConstPtr mask, double -GetMapMedian(G3SkyMapConstPtr m, G3SkyMapConstPtr mask, bool ignore_zeros, bool ignore_nans, +GetMapMedian(G3SkyMapConstPtr m, G3SkyMapMaskConstPtr 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()); + size_t npix = mask ? mask->sum() : (ignore_zeros ? m->NpixAllocated() : m->size()); if (npix == 0) return 0; @@ -473,7 +439,7 @@ GetMapMedian(G3SkyMapConstPtr m, G3SkyMapConstPtr mask, bool ignore_zeros, bool std::vector -GetMapMinMax(G3SkyMapConstPtr m, G3SkyMapConstPtr mask, bool ignore_zeros, bool ignore_nans, +GetMapMinMax(G3SkyMapConstPtr m, G3SkyMapMaskConstPtr mask, bool ignore_zeros, bool ignore_nans, bool ignore_infs) { double min = 0.0 / 0.0; @@ -500,7 +466,7 @@ GetMapMinMax(G3SkyMapConstPtr m, G3SkyMapConstPtr mask, bool ignore_zeros, bool 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 +578,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 " @@ -694,7 +649,7 @@ void maputils_pybindings(void){ "copied from the input map."); bp::def("get_map_stats", GetMapStats, - (bp::arg("map"), bp::arg("mask")=G3SkyMapConstPtr(), bp::arg("order")=2, + (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 " @@ -703,21 +658,21 @@ void maputils_pybindings(void){ "the non-zero pixels in the mask are included."); bp::def("get_map_median", GetMapMedian, - (bp::arg("map"), bp::arg("mask")=G3SkyMapConstPtr(), + (bp::arg("map"), bp::arg("mask")=G3SkyMapMaskConstPtr(), 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("map"), bp::arg("mask")=G3SkyMapMaskConstPtr(), 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 +684,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..9c8eade9 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(numpy.where(numpy.asarray(mask.to_map()).ravel())[0]) ^ masked) diff --git a/maps/tests/healpix_maps.py b/maps/tests/healpix_maps.py index 3d228785..3e926bfe 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(np.where(np.asarray(mask.to_map()).ravel())[0]) ^ masked) x.shift_ra = True 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) From 70afcd7f58993415c5a60047118cce643766ef63 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Thu, 14 Jul 2022 11:53:11 -0500 Subject: [PATCH 25/55] apply_mask method for G3SkyMaps --- maps/include/maps/G3SkyMap.h | 2 ++ maps/src/G3SkyMap.cxx | 18 ++++++++++++++++++ 2 files changed, 20 insertions(+) diff --git a/maps/include/maps/G3SkyMap.h b/maps/include/maps/G3SkyMap.h index 958cfe35..9127090e 100644 --- a/maps/include/maps/G3SkyMap.h +++ b/maps/include/maps/G3SkyMap.h @@ -127,6 +127,8 @@ class G3SkyMap { virtual bool all() const; virtual bool any() const; + virtual void ApplyMask(const G3SkyMapMask &mask, bool inverse=false); + // Pointing information std::vector AnglesToPixels(const std::vector & alphas, const std::vector & deltas) const; diff --git a/maps/src/G3SkyMap.cxx b/maps/src/G3SkyMap.cxx index c7ae6de4..a4f1064b 100644 --- a/maps/src/G3SkyMap.cxx +++ b/maps/src/G3SkyMap.cxx @@ -693,6 +693,19 @@ pyskymap_bool(G3SkyMap &skymap) return false; } +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; + } +} + G3SkyMapWeights & G3SkyMapWeights::operator+=(const G3SkyMapWeights &rhs) { @@ -1104,6 +1117,11 @@ PYBINDINGS("maps") { "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.") + .def(bp::self += bp::self) .def(bp::self *= bp::self) .def(bp::self -= bp::self) From 3b75d39a0f9202c32e5e1e2a2c4d80aa856179a0 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Wed, 20 Jul 2022 02:05:36 -0500 Subject: [PATCH 26/55] apply_mask for weights consistent handling of polarized weight attributes --- maps/include/maps/G3SkyMap.h | 2 + maps/python/skymapaddons.py | 9 +++- maps/src/G3SkyMap.cxx | 84 ++++++++++++++++++++++++------------ 3 files changed, 66 insertions(+), 29 deletions(-) diff --git a/maps/include/maps/G3SkyMap.h b/maps/include/maps/G3SkyMap.h index 9127090e..c29519e1 100644 --- a/maps/include/maps/G3SkyMap.h +++ b/maps/include/maps/G3SkyMap.h @@ -399,6 +399,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/python/skymapaddons.py b/maps/python/skymapaddons.py index 2b595fdf..8055672e 100644 --- a/maps/python/skymapaddons.py +++ b/maps/python/skymapaddons.py @@ -160,12 +160,17 @@ 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 diff --git a/maps/src/G3SkyMap.cxx b/maps/src/G3SkyMap.cxx index a4f1064b..dddbca1a 100644 --- a/maps/src/G3SkyMap.cxx +++ b/maps/src/G3SkyMap.cxx @@ -706,20 +706,40 @@ G3SkyMap::ApplyMask(const G3SkyMapMask &mask, bool inverse) } } +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; } @@ -727,14 +747,18 @@ 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; \ } @@ -914,13 +938,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(); @@ -933,20 +960,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; } @@ -955,14 +974,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") { @@ -1120,7 +1143,8 @@ PYBINDINGS("maps") { .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.") + "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) @@ -1208,6 +1232,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.") From 6889f6826b516cb68f455a6a8fde34e7fded1e91 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Wed, 20 Jul 2022 11:07:09 -0500 Subject: [PATCH 27/55] rename get_map_stats -> get_map_moments --- maps/include/maps/maputils.h | 2 +- maps/src/maputils.cxx | 4 ++-- maps/tests/flatskymap_operators.py | 8 ++++---- 3 files changed, 7 insertions(+), 7 deletions(-) diff --git a/maps/include/maps/maputils.h b/maps/include/maps/maputils.h index 8ce46e85..2510de13 100644 --- a/maps/include/maps/maputils.h +++ b/maps/include/maps/maputils.h @@ -36,7 +36,7 @@ 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, G3SkyMapMaskConstPtr mask=NULL, int order=2, +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 median of the input map, optionally excluding any pixels that are diff --git a/maps/src/maputils.cxx b/maps/src/maputils.cxx index 620c2ad9..703a698f 100644 --- a/maps/src/maputils.cxx +++ b/maps/src/maputils.cxx @@ -348,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, G3SkyMapMaskConstPtr mask, +std::vector GetMapMoments(G3SkyMapConstPtr m, G3SkyMapMaskConstPtr mask, int order, bool ignore_zeros, bool ignore_nans, bool ignore_infs) { size_t n = 0; @@ -648,7 +648,7 @@ 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::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 " diff --git a/maps/tests/flatskymap_operators.py b/maps/tests/flatskymap_operators.py index d92b4f3d..63a5bb8d 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, get_map_median from scipy.stats import skew, kurtosis # Sparse extension operators @@ -253,19 +253,19 @@ # 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) 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) assert(np.allclose(med2, med0)) np.asarray(mpad)[np.asarray(mpad) == 0] = np.nan - stats3 = get_map_stats(mpad, order=4, ignore_nans=True) + stats3 = get_map_moments(mpad, order=4, ignore_nans=True) assert(np.allclose(stats3, stats0)) med3 = get_map_median(mpad, ignore_nans=True) assert(np.allclose(med3, med0)) From 982e3ec59a5648203edd18bf29a439ae85b001c0 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Sat, 30 Jul 2022 16:22:53 -0500 Subject: [PATCH 28/55] clone, element-wise comparison operators for masks --- maps/include/maps/G3SkyMapMask.h | 5 +++++ maps/src/G3SkyMapMask.cxx | 38 ++++++++++++++++++++++++++++++++ 2 files changed, 43 insertions(+) diff --git a/maps/include/maps/G3SkyMapMask.h b/maps/include/maps/G3SkyMapMask.h index b6564383..222b28b3 100644 --- a/maps/include/maps/G3SkyMapMask.h +++ b/maps/include/maps/G3SkyMapMask.h @@ -26,6 +26,9 @@ class G3SkyMapMask : public G3FrameObject { G3SkyMapMask(const G3SkyMapMask &); virtual ~G3SkyMapMask() {}; + // Copy + boost::shared_ptr Clone(bool copy_data = true) const; + // Return a (modifiable) pixel value std::vector::reference operator [] (size_t i); bool operator [] (size_t i) const { return this->at(i); }; @@ -46,6 +49,8 @@ class G3SkyMapMask : public G3FrameObject { 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; // Information bool IsCompatible(const G3SkyMap &map) const; diff --git a/maps/src/G3SkyMapMask.cxx b/maps/src/G3SkyMapMask.cxx index a798f298..ad3189ce 100644 --- a/maps/src/G3SkyMapMask.cxx +++ b/maps/src/G3SkyMapMask.cxx @@ -30,6 +30,15 @@ G3SkyMapMask::G3SkyMapMask(const G3SkyMapMask &m) : G3FrameObject() 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()); +} + std::vector::reference G3SkyMapMask::operator [] (size_t i) { @@ -170,6 +179,30 @@ G3SkyMapMask::operator ^(const G3SkyMapMask &rhs) const return mask; } +G3SkyMapMask +G3SkyMapMask::operator ==(const G3SkyMapMask &rhs) const +{ + g3_assert(IsCompatible(rhs)); + + G3SkyMapMask mask(*Parent()); + for (size_t i = 0; i < data_.size(); i++) + mask[i] = at(i) == rhs.at(i); + + return mask; +} + +G3SkyMapMask +G3SkyMapMask::operator !=(const G3SkyMapMask &rhs) const +{ + g3_assert(IsCompatible(rhs)); + + G3SkyMapMask mask(*Parent()); + for (size_t i = 0; i < data_.size(); i++) + mask[i] = at(i) != rhs.at(i); + + return mask; +} + bool G3SkyMapMask::IsCompatible(const G3SkyMap &map) const { @@ -323,6 +356,9 @@ PYBINDINGS("maps") "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("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_readonly("parent", &G3SkyMapMask::Parent, "\"Parent\" map which " "contains no data, but can be used to retrieve the parameters of " "the map to which this mask corresponds.") @@ -343,6 +379,8 @@ PYBINDINGS("maps") .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(); From b5faccaa99fbb5657691a3cf29b4b30947436257 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Sat, 30 Jul 2022 16:24:24 -0500 Subject: [PATCH 29/55] get parent attributes from mask --- maps/python/skymapaddons.py | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/maps/python/skymapaddons.py b/maps/python/skymapaddons.py index 8055672e..125720c4 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 @@ -152,7 +152,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__ @@ -175,6 +177,12 @@ def skymapweights_setattr(self, 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(): From ad343e8bdbd98482cbf4a3979815379a3f9550c3 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Sat, 30 Jul 2022 16:24:39 -0500 Subject: [PATCH 30/55] mask tests --- maps/CMakeLists.txt | 1 + maps/tests/mask_operators.py | 73 ++++++++++++++++++++++++++++++++++++ 2 files changed, 74 insertions(+) create mode 100755 maps/tests/mask_operators.py diff --git a/maps/CMakeLists.txt b/maps/CMakeLists.txt index 8ddf08e5..1edd89d9 100644 --- a/maps/CMakeLists.txt +++ b/maps/CMakeLists.txt @@ -28,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/tests/mask_operators.py b/maps/tests/mask_operators.py new file mode 100755 index 00000000..279878dd --- /dev/null +++ b/maps/tests/mask_operators.py @@ -0,0 +1,73 @@ +#!/usr/bin/env python + +import numpy as np +from spt3g import core, maps + +shape = (300, 300) +x = maps.FlatSkyMap(np.random.randn(*shape), core.G3Units.arcmin) + +# float comparison +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 +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 +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 +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_pixels() +xpix, _ = x4.nonzero_pixels() + +assert(len(set(mpix) ^ set(xpix)) == 0) From 9c0ef00fa94952cdf729d5666d27759c6ef2f106 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Sat, 30 Jul 2022 16:47:39 -0500 Subject: [PATCH 31/55] name unused --- maps/src/G3SkyMap.cxx | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/maps/src/G3SkyMap.cxx b/maps/src/G3SkyMap.cxx index dddbca1a..638c3d4d 100644 --- a/maps/src/G3SkyMap.cxx +++ b/maps/src/G3SkyMap.cxx @@ -625,7 +625,7 @@ pyskymap_pow(const G3SkyMap &a, const G3SkyMap &b) return rv; } -#define skymap_comp(name, oper) \ +#define skymap_comp(oper) \ G3SkyMapMask \ G3SkyMap::operator oper(const G3SkyMap &rhs) \ { \ @@ -639,14 +639,14 @@ G3SkyMap::operator oper(const G3SkyMap &rhs) \ 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) \ +#define skymap_compd(oper) \ G3SkyMapMask \ G3SkyMap::operator oper(double rhs) \ { \ @@ -658,12 +658,12 @@ G3SkyMap::operator oper(double rhs) \ 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(>) bool G3SkyMap::all() const From 4d10a753e1f7abc358a3a1d15e01aac007c0cc52 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Sat, 30 Jul 2022 18:05:52 -0500 Subject: [PATCH 32/55] iterator for masks --- maps/include/maps/G3SkyMapMask.h | 37 ++++++++++++ maps/src/G3SkyMap.cxx | 7 ++- maps/src/G3SkyMapMask.cxx | 100 ++++++++++++++++++++----------- 3 files changed, 106 insertions(+), 38 deletions(-) diff --git a/maps/include/maps/G3SkyMapMask.h b/maps/include/maps/G3SkyMapMask.h index 222b28b3..6694da35 100644 --- a/maps/include/maps/G3SkyMapMask.h +++ b/maps/include/maps/G3SkyMapMask.h @@ -43,6 +43,7 @@ class G3SkyMapMask : public G3FrameObject { bool all() const; bool any() const; size_t sum() const; + size_t size() const; std::vector NonZeroPixels() const; G3SkyMapMask operator ~() const; @@ -63,6 +64,42 @@ class G3SkyMapMask : public G3FrameObject { // 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); diff --git a/maps/src/G3SkyMap.cxx b/maps/src/G3SkyMap.cxx index 638c3d4d..93a7e4d5 100644 --- a/maps/src/G3SkyMap.cxx +++ b/maps/src/G3SkyMap.cxx @@ -532,10 +532,11 @@ 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 (size_t i = 0; i < a.size(); i++) { - if (b.at(i) && a.at(i) != 0) - (*rv)[i] = a.at(i); + for (auto i: b) { + if (b.at(i.first) && a.at(i.first) != 0) + (*rv)[i.first] = a.at(i.first); } return rv; diff --git a/maps/src/G3SkyMapMask.cxx b/maps/src/G3SkyMapMask.cxx index ad3189ce..975e5a52 100644 --- a/maps/src/G3SkyMapMask.cxx +++ b/maps/src/G3SkyMapMask.cxx @@ -39,6 +39,28 @@ G3SkyMapMask::Clone(bool copy_data) const return boost::make_shared(*Parent()); } +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) { @@ -48,6 +70,8 @@ G3SkyMapMask::operator [] (size_t i) bool G3SkyMapMask::at (size_t i) const { + if (i >= data_.size()) + return false; return data_[i]; } @@ -56,8 +80,8 @@ G3SkyMapMask::operator |=(const G3SkyMapMask &rhs) { g3_assert(IsCompatible(rhs)); - for (size_t i = 0; i < data_.size(); i++) - (*this)[i] = rhs.at(i) || (*this)[i]; + for (size_t i = 0; i < size(); i++) + (*this)[i] = rhs.at(i) || at(i); return *this; } @@ -67,8 +91,8 @@ G3SkyMapMask::operator &=(const G3SkyMapMask &rhs) { g3_assert(IsCompatible(rhs)); - for (size_t i = 0; i < data_.size(); i++) - (*this)[i] = rhs.at(i) && (*this)[i]; + for (auto i: *this) + (*this)[i.first] = rhs.at(i.first) && i.second; return *this; } @@ -78,8 +102,8 @@ G3SkyMapMask::operator ^=(const G3SkyMapMask &rhs) { g3_assert(IsCompatible(rhs)); - for (size_t i = 0; i < data_.size(); i++) - (*this)[i] = rhs.at(i) ^ (*this)[i]; + for (size_t i = 0; i < size(); i++) + (*this)[i] = rhs.at(i) ^ at(i); return *this; } @@ -87,8 +111,8 @@ G3SkyMapMask::operator ^=(const G3SkyMapMask &rhs) G3SkyMapMask & G3SkyMapMask::invert() { - for (size_t i = 0; i < data_.size(); i++) - (*this)[i] = !(*this)[i]; + for (size_t i = 0; i < size(); i++) + (*this)[i] = !at(i); return *this; } @@ -96,7 +120,7 @@ G3SkyMapMask::invert() bool G3SkyMapMask::all() const { - for (size_t i = 0; i < data_.size(); i++) + for (size_t i = 0; i < size(); i++) if (at(i) == 0) return false; return true; @@ -105,8 +129,8 @@ G3SkyMapMask::all() const bool G3SkyMapMask::any() const { - for (size_t i = 0; i < data_.size(); i++) - if (at(i) != 0) + for (auto i: *this) + if (i.second != 0) return true; return false; } @@ -115,8 +139,8 @@ size_t G3SkyMapMask::sum() const { size_t s = 0; - for (size_t i = 0; i < data_.size(); i++) - s += at(i); + for (auto i: *this) + s += i.second; return s; } @@ -125,9 +149,9 @@ G3SkyMapMask::NonZeroPixels() const { std::vector indices; - for (size_t i = 0; i < data_.size(); i++) { - if (at(i)) - indices.push_back(i); + for (auto i: *this) { + if (i.second) + indices.push_back(i.first); } return indices; @@ -137,8 +161,9 @@ G3SkyMapMask G3SkyMapMask::operator ~() const { G3SkyMapMask mask(*Parent()); - for (size_t i = 0; i < data_.size(); i++) - mask[i] = !at(i); + for (size_t i = 0; i < size(); i++) + if (!at(i)) + mask[i] = true; return mask; } @@ -149,8 +174,9 @@ G3SkyMapMask::operator |(const G3SkyMapMask &rhs) const g3_assert(IsCompatible(rhs)); G3SkyMapMask mask(*Parent()); - for (size_t i = 0; i < data_.size(); i++) - mask[i] = at(i) || rhs.at(i); + for (size_t i = 0; i < size(); i++) + if (at(i) || rhs.at(i)) + mask[i] = true; return mask; } @@ -161,8 +187,9 @@ G3SkyMapMask::operator &(const G3SkyMapMask &rhs) const g3_assert(IsCompatible(rhs)); G3SkyMapMask mask(*Parent()); - for (size_t i = 0; i < data_.size(); i++) - mask[i] = at(i) && rhs.at(i); + for (auto i: *this) + if (i.second && rhs.at(i.first)) + mask[i.first] = true; return mask; } @@ -173,8 +200,9 @@ G3SkyMapMask::operator ^(const G3SkyMapMask &rhs) const g3_assert(IsCompatible(rhs)); G3SkyMapMask mask(*Parent()); - for (size_t i = 0; i < data_.size(); i++) - mask[i] = at(i) ^ rhs.at(i); + for (size_t i = 0; i < size(); i++) + if (at(i) ^ rhs.at(i)) + mask[i] = true; return mask; } @@ -185,8 +213,9 @@ G3SkyMapMask::operator ==(const G3SkyMapMask &rhs) const g3_assert(IsCompatible(rhs)); G3SkyMapMask mask(*Parent()); - for (size_t i = 0; i < data_.size(); i++) - mask[i] = at(i) == rhs.at(i); + for (size_t i = 0; i < size(); i++) + if (at(i) == rhs.at(i)) + mask[i] = true; return mask; } @@ -197,8 +226,9 @@ G3SkyMapMask::operator !=(const G3SkyMapMask &rhs) const g3_assert(IsCompatible(rhs)); G3SkyMapMask mask(*Parent()); - for (size_t i = 0; i < data_.size(); i++) - mask[i] = at(i) != rhs.at(i); + for (size_t i = 0; i < size(); i++) + if (at(i) != rhs.at(i)) + mask[i] = true; return mask; } @@ -220,9 +250,9 @@ G3SkyMapMask::MakeBinaryMap() const { G3SkyMapPtr out = Parent()->Clone(); - for (size_t i = 0; i < data_.size(); i++) { - if (at(i)) - (*out)[i] = 1.0; + for (auto i: *this) { + if (i.second) + (*out)[i.first] = 1.0; } return out; @@ -249,7 +279,7 @@ skymapmask_getitem(const G3SkyMapMask &m, boost::python::object index) i = extract(index)(); if (i < 0) - i = m.Parent()->size() + i; + i = m.size() + i; } else if (extract(index).check()) { int x, y; @@ -280,7 +310,7 @@ skymapmask_getitem(const G3SkyMapMask &m, boost::python::object index) boost::python::throw_error_already_set(); } - if (i < 0 || size_t(i) >= m.Parent()->size()) { + if (i < 0 || size_t(i) >= m.size()) { PyErr_SetString(PyExc_IndexError, "Index out of range"); boost::python::throw_error_already_set(); } @@ -298,7 +328,7 @@ skymapmask_setitem(G3SkyMapMask &m, boost::python::object index, bool val) i = extract(index)(); if (i < 0) - i = m.Parent()->size() + i; + i = m.size() + i; } else if (extract(index).check()) { int x, y; @@ -329,7 +359,7 @@ skymapmask_setitem(G3SkyMapMask &m, boost::python::object index, bool val) boost::python::throw_error_already_set(); } - if (i < 0 || size_t(i) >= m.Parent()->size()) { + if (i < 0 || size_t(i) >= m.size()) { PyErr_SetString(PyExc_IndexError, "Index out of range"); boost::python::throw_error_already_set(); } From 1d60036d71977bd43b8f1a0e8bb8292b564da71c Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Sat, 30 Jul 2022 19:39:49 -0500 Subject: [PATCH 33/55] efficient apply_mask functions --- maps/include/maps/FlatSkyMap.h | 1 + maps/include/maps/HealpixSkyMap.h | 2 ++ maps/src/FlatSkyMap.cxx | 8 ++++++++ maps/src/HealpixSkyMap.cxx | 8 ++++++++ 4 files changed, 19 insertions(+) diff --git a/maps/include/maps/FlatSkyMap.h b/maps/include/maps/FlatSkyMap.h index 9a5f537b..7f97e50f 100644 --- a/maps/include/maps/FlatSkyMap.h +++ b/maps/include/maps/FlatSkyMap.h @@ -93,6 +93,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); diff --git a/maps/include/maps/HealpixSkyMap.h b/maps/include/maps/HealpixSkyMap.h index e07a88af..b2613bc6 100644 --- a/maps/include/maps/HealpixSkyMap.h +++ b/maps/include/maps/HealpixSkyMap.h @@ -72,6 +72,8 @@ class HealpixSkyMap : public G3FrameObject, public G3SkyMap { 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();} double res() const; diff --git a/maps/src/FlatSkyMap.cxx b/maps/src/FlatSkyMap.cxx index 58ee8710..81ebfdba 100644 --- a/maps/src/FlatSkyMap.cxx +++ b/maps/src/FlatSkyMap.cxx @@ -636,6 +636,14 @@ FlatSkyMap::NonZeroPixels(std::vector &indices, } } +void +FlatSkyMap::ApplyMask(const G3SkyMapMask &mask, bool inverse) +{ + 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_}; } diff --git a/maps/src/HealpixSkyMap.cxx b/maps/src/HealpixSkyMap.cxx index bb756ef3..0682bea7 100644 --- a/maps/src/HealpixSkyMap.cxx +++ b/maps/src/HealpixSkyMap.cxx @@ -859,6 +859,14 @@ HealpixSkyMap::NonZeroPixels(std::vector &indices, } } +void +HealpixSkyMap::ApplyMask(const G3SkyMapMask &mask, bool inverse) +{ + for (auto i: *this) + if (i.second != 0 && mask.at(i.first) == inverse) + (*this)[i.first] = 0; +} + std::vector HealpixSkyMap::shape() const From 4f97e714eaf1fa2d4a09d8ada86696c46d5092b7 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Sat, 30 Jul 2022 19:46:23 -0500 Subject: [PATCH 34/55] efficient any() functions --- maps/include/maps/FlatSkyMap.h | 1 + maps/include/maps/HealpixSkyMap.h | 2 +- maps/src/FlatSkyMap.cxx | 8 ++++++++ maps/src/HealpixSkyMap.cxx | 8 ++++++++ 4 files changed, 18 insertions(+), 1 deletion(-) diff --git a/maps/include/maps/FlatSkyMap.h b/maps/include/maps/FlatSkyMap.h index 7f97e50f..793c049f 100644 --- a/maps/include/maps/FlatSkyMap.h +++ b/maps/include/maps/FlatSkyMap.h @@ -94,6 +94,7 @@ class FlatSkyMap : public G3FrameObject, public G3SkyMap { void NonZeroPixels(std::vector &indices, std::vector &data) const; void ApplyMask(const G3SkyMapMask &mask, bool inverse=false) override; + bool any() const override; void SetProj(MapProjection proj); void SetAlphaCenter(double alpha); diff --git a/maps/include/maps/HealpixSkyMap.h b/maps/include/maps/HealpixSkyMap.h index b2613bc6..fc75c359 100644 --- a/maps/include/maps/HealpixSkyMap.h +++ b/maps/include/maps/HealpixSkyMap.h @@ -71,8 +71,8 @@ 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; + bool any() const override; size_t nside() const {return info_.nside();} bool nested() const {return info_.nested();} diff --git a/maps/src/FlatSkyMap.cxx b/maps/src/FlatSkyMap.cxx index 81ebfdba..f852893a 100644 --- a/maps/src/FlatSkyMap.cxx +++ b/maps/src/FlatSkyMap.cxx @@ -664,6 +664,14 @@ size_t FlatSkyMap::NpixNonZero() const { return 0; } +bool FlatSkyMap::any() const +{ + for (auto i: *this) + if (i.second != 0) + return true; + return false; +} + #define GETSET(name, cname, type) \ type FlatSkyMap::name() const \ { \ diff --git a/maps/src/HealpixSkyMap.cxx b/maps/src/HealpixSkyMap.cxx index 0682bea7..1f2999d8 100644 --- a/maps/src/HealpixSkyMap.cxx +++ b/maps/src/HealpixSkyMap.cxx @@ -910,6 +910,14 @@ size_t HealpixSkyMap::NpixNonZero() const return sz; } +bool HealpixSkyMap::any() const +{ + for (auto i: *this) + if (i.second != 0) + return true; + return false; +} + size_t HealpixSkyMap::AngleToPixel(double alpha, double delta) const { From a73163529cfa5d69a7dac79d5404f545385ab6c8 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Sat, 30 Jul 2022 20:05:16 -0500 Subject: [PATCH 35/55] mask assignment test --- maps/tests/mask_operators.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/maps/tests/mask_operators.py b/maps/tests/mask_operators.py index 279878dd..16a06b93 100755 --- a/maps/tests/mask_operators.py +++ b/maps/tests/mask_operators.py @@ -6,6 +6,14 @@ shape = (300, 300) x = maps.FlatSkyMap(np.random.randn(*shape), core.G3Units.arcmin) +# element-wise assignment +m = x.to_mask() +v = m[5, 5] +m[5, 5] = not v +assert(m[5, 5] != v) +m[5, 5] = v +assert(m[5, 5] == v) + # float comparison m1 = x == 0 m2 = x != 0 From c8ccba01cb372c3822b3f081ccccf089a4c64b8b Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Sat, 30 Jul 2022 21:08:21 -0500 Subject: [PATCH 36/55] test parent attributes --- maps/src/G3SkyMapMask.cxx | 3 ++- maps/tests/mask_operators.py | 8 ++++++-- 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/maps/src/G3SkyMapMask.cxx b/maps/src/G3SkyMapMask.cxx index 975e5a52..5be48170 100644 --- a/maps/src/G3SkyMapMask.cxx +++ b/maps/src/G3SkyMapMask.cxx @@ -389,9 +389,10 @@ PYBINDINGS("maps") .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_readonly("parent", &G3SkyMapMask::Parent, "\"Parent\" map which " + .add_property("parent", &G3SkyMapMask::Parent, "\"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) diff --git a/maps/tests/mask_operators.py b/maps/tests/mask_operators.py index 16a06b93..4b2dcc23 100755 --- a/maps/tests/mask_operators.py +++ b/maps/tests/mask_operators.py @@ -3,11 +3,15 @@ import numpy as np from spt3g import core, maps -shape = (300, 300) +shape = (300, 500) x = maps.FlatSkyMap(np.random.randn(*shape), core.G3Units.arcmin) +m = x.to_mask() + +# attributes +assert(m.size == x.size) +assert(m.shape == x.shape) # element-wise assignment -m = x.to_mask() v = m[5, 5] m[5, 5] = not v assert(m[5, 5] != v) From 8cb1975fffc68c8fd1dc77d52e82c03c9b0c0197 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Sat, 30 Jul 2022 21:54:08 -0500 Subject: [PATCH 37/55] use mask class in MapTODMasker --- maps/src/MapTODMasker.cxx | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/maps/src/MapTODMasker.cxx b/maps/src/MapTODMasker.cxx index 8fbc0e6f..d94036a3 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) @@ -105,12 +105,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); + mask_->Parent()->coord_ref, alpha, delta); + auto detpointing = mask_->Parent()->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; From 7c88153fd6f3d517de331a27e81b5f94085c739e Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Sat, 30 Jul 2022 23:40:44 -0500 Subject: [PATCH 38/55] apply inverse mask in-place in masks too --- maps/include/maps/G3SkyMapMask.h | 2 ++ maps/src/G3SkyMapMask.cxx | 15 +++++++++++++++ maps/tests/mask_operators.py | 3 +++ 3 files changed, 20 insertions(+) diff --git a/maps/include/maps/G3SkyMapMask.h b/maps/include/maps/G3SkyMapMask.h index 6694da35..cf222272 100644 --- a/maps/include/maps/G3SkyMapMask.h +++ b/maps/include/maps/G3SkyMapMask.h @@ -53,6 +53,8 @@ class G3SkyMapMask : public G3FrameObject { 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; diff --git a/maps/src/G3SkyMapMask.cxx b/maps/src/G3SkyMapMask.cxx index 5be48170..4ea9f04e 100644 --- a/maps/src/G3SkyMapMask.cxx +++ b/maps/src/G3SkyMapMask.cxx @@ -233,6 +233,16 @@ G3SkyMapMask::operator !=(const G3SkyMapMask &rhs) const 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 { @@ -403,6 +413,11 @@ PYBINDINGS("maps") .def("sum", &G3SkyMapMask::sum, "Sum of all elements in mask") .def("nonzero_pixels", &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) diff --git a/maps/tests/mask_operators.py b/maps/tests/mask_operators.py index 4b2dcc23..2b45aeb5 100755 --- a/maps/tests/mask_operators.py +++ b/maps/tests/mask_operators.py @@ -83,3 +83,6 @@ xpix, _ = x4.nonzero_pixels() assert(len(set(mpix) ^ set(xpix)) == 0) + +m1.apply_mask(m2) +assert(not m1.any()) From 7fc7b5663d63150e2674454e73b832464f6229a5 Mon Sep 17 00:00:00 2001 From: Nathan Whitehorn Date: Mon, 1 Aug 2022 11:15:32 -0400 Subject: [PATCH 39/55] Fix test failure. The fix isn't great, but c'est la vie. --- maps/src/G3SkyMapMask.cxx | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/maps/src/G3SkyMapMask.cxx b/maps/src/G3SkyMapMask.cxx index 4ea9f04e..8e7bfb64 100644 --- a/maps/src/G3SkyMapMask.cxx +++ b/maps/src/G3SkyMapMask.cxx @@ -386,6 +386,14 @@ skymapmask_pyinvert(G3SkyMapMaskPtr m) 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()); +} + PYBINDINGS("maps") { using namespace boost::python; @@ -399,7 +407,7 @@ PYBINDINGS("maps") .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.") - .add_property("parent", &G3SkyMapMask::Parent, "\"Parent\" map which " + .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") From a06dc9cfe9ede12c36a2c0a74286a5a483bdeb22 Mon Sep 17 00:00:00 2001 From: Nathan Whitehorn Date: Mon, 1 Aug 2022 11:16:02 -0400 Subject: [PATCH 40/55] Regularize implicit pointer conversions. --- maps/src/FlatSkyMap.cxx | 2 -- maps/src/G3SkyMap.cxx | 1 + maps/src/HealpixSkyMap.cxx | 2 -- 3 files changed, 1 insertion(+), 4 deletions(-) diff --git a/maps/src/FlatSkyMap.cxx b/maps/src/FlatSkyMap.cxx index f852893a..112f1c15 100644 --- a/maps/src/FlatSkyMap.cxx +++ b/maps/src/FlatSkyMap.cxx @@ -1253,8 +1253,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 93a7e4d5..70f0847d 100644 --- a/maps/src/G3SkyMap.cxx +++ b/maps/src/G3SkyMap.cxx @@ -1196,6 +1196,7 @@ PYBINDINGS("maps") { .def("any", &G3SkyMap::any) .def("all", &G3SkyMap::all) ; + boost::python::implicitly_convertible(); EXPORT_FRAMEOBJECT(G3SkyMapWeights, init<>(), "Polarized (Mueller matrix) or unpolarized (scalar) map pixel weights.") diff --git a/maps/src/HealpixSkyMap.cxx b/maps/src/HealpixSkyMap.cxx index 1f2999d8..23003d61 100644 --- a/maps/src/HealpixSkyMap.cxx +++ b/maps/src/HealpixSkyMap.cxx @@ -1202,8 +1202,6 @@ PYBINDINGS("maps") #endif implicitly_convertible(); - implicitly_convertible(); implicitly_convertible(); - implicitly_convertible(); } From caa57afc5b40ab12265e8138f911fdfc66caf6ef Mon Sep 17 00:00:00 2001 From: Nathan Whitehorn Date: Mon, 1 Aug 2022 11:49:02 -0400 Subject: [PATCH 41/55] Ask masked get and set operations for FlatSkyMap. --- maps/src/FlatSkyMap.cxx | 41 ++++++++++++++++++++++++++++++++++++ maps/tests/mask_operators.py | 13 ++++++++++++ 2 files changed, 54 insertions(+) diff --git a/maps/src/FlatSkyMap.cxx b/maps/src/FlatSkyMap.cxx index 112f1c15..b1bce27c 100644 --- a/maps/src/FlatSkyMap.cxx +++ b/maps/src/FlatSkyMap.cxx @@ -1044,6 +1044,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) { @@ -1241,6 +1280,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(); diff --git a/maps/tests/mask_operators.py b/maps/tests/mask_operators.py index 2b45aeb5..ea81fb1a 100755 --- a/maps/tests/mask_operators.py +++ b/maps/tests/mask_operators.py @@ -18,6 +18,19 @@ m[5, 5] = v assert(m[5, 5] == v) +# masked assignment and retrieval +pospixels = x[x > 0] +assert(len(pospixels) == (x > 0).sum()) +assert((np.asarray(pospixels) > 0).all()) + +x2 = maps.FlatSkyMap(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) + # float comparison m1 = x == 0 m2 = x != 0 From 025a9021961a132218d89831ae53235769282cf2 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Mon, 1 Aug 2022 12:59:14 -0500 Subject: [PATCH 42/55] Add masked get and set operations for HealpixSkyMap, test all the things --- maps/src/HealpixSkyMap.cxx | 42 ++++++++ maps/tests/mask_operators.py | 203 ++++++++++++++++++----------------- 2 files changed, 149 insertions(+), 96 deletions(-) diff --git a/maps/src/HealpixSkyMap.cxx b/maps/src/HealpixSkyMap.cxx index 23003d61..0b7d9170 100644 --- a/maps/src/HealpixSkyMap.cxx +++ b/maps/src/HealpixSkyMap.cxx @@ -1069,6 +1069,45 @@ HealpixSkyMap_nonzeropixels(const HealpixSkyMap &m) return boost::python::make_tuple(i, d); } +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) { @@ -1187,6 +1226,9 @@ PYBINDINGS("maps") "holes or very small filling factors. " "If set to True, converts the map to this representation." ) + .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.") diff --git a/maps/tests/mask_operators.py b/maps/tests/mask_operators.py index ea81fb1a..dd700158 100755 --- a/maps/tests/mask_operators.py +++ b/maps/tests/mask_operators.py @@ -3,99 +3,110 @@ import numpy as np from spt3g import core, maps -shape = (300, 500) -x = maps.FlatSkyMap(np.random.randn(*shape), core.G3Units.arcmin) -m = x.to_mask() - -# attributes -assert(m.size == x.size) -assert(m.shape == x.shape) - -# element-wise assignment -v = m[5, 5] -m[5, 5] = not v -assert(m[5, 5] != v) -m[5, 5] = v -assert(m[5, 5] == v) - -# masked assignment and retrieval -pospixels = x[x > 0] -assert(len(pospixels) == (x > 0).sum()) -assert((np.asarray(pospixels) > 0).all()) - -x2 = maps.FlatSkyMap(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) - -# float comparison -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 -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 -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 -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_pixels() -xpix, _ = x4.nonzero_pixels() - -assert(len(set(mpix) ^ set(xpix)) == 0) - -m1.apply_mask(m2) -assert(not m1.any()) +maplist = [ + maps.FlatSkyMap(np.random.randn(300, 500), core.G3Units.arcmin), + maps.HealpixSkyMap(np.random.randn(12 * 256 * 256)), +] + +for x in maplist: + 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 + 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) + + # float comparison + 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 + 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 + 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 + 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_pixels() + xpix, _ = x4.nonzero_pixels() + + assert(len(set(mpix) ^ set(xpix)) == 0) + + m1.apply_mask(m2) + assert(not m1.any()) From 70907885b21153d8ab2fd7cd90944c3d37be22dc Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Mon, 1 Aug 2022 15:28:06 -0500 Subject: [PATCH 43/55] need to add base class get/set methods back in --- maps/src/HealpixSkyMap.cxx | 30 ++++++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/maps/src/HealpixSkyMap.cxx b/maps/src/HealpixSkyMap.cxx index 0b7d9170..f41c350d 100644 --- a/maps/src/HealpixSkyMap.cxx +++ b/maps/src/HealpixSkyMap.cxx @@ -1069,6 +1069,34 @@ 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) { @@ -1226,6 +1254,8 @@ 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) From 820226085a0793f6de2d08047908408bae0e2586 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Mon, 1 Aug 2022 16:17:34 -0500 Subject: [PATCH 44/55] masks have no units --- maps/src/G3SkyMapMask.cxx | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/maps/src/G3SkyMapMask.cxx b/maps/src/G3SkyMapMask.cxx index 8e7bfb64..e7d6eca5 100644 --- a/maps/src/G3SkyMapMask.cxx +++ b/maps/src/G3SkyMapMask.cxx @@ -10,7 +10,13 @@ G3SkyMapMask::G3SkyMapMask(const G3SkyMap &parent, bool use_data, bool zero_nans, bool zero_infs) : G3FrameObject() { - parent_ = parent.Clone(false); + // 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) { for (size_t i = 0; i < parent.size(); i++) { From 5a38ebd4a524d5975d7cac38633033b11b10dc5c Mon Sep 17 00:00:00 2001 From: Alexandra Rahlin Date: Tue, 2 Aug 2022 18:08:47 -0500 Subject: [PATCH 45/55] Add mask documentation to maps readme --- maps/README.rst | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/maps/README.rst b/maps/README.rst index 0c4ae86c..155b68bd 100644 --- a/maps/README.rst +++ b/maps/README.rst @@ -58,6 +58,24 @@ By default, both Healpix and flat-sky maps are initialized in sparse mode. This 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. +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. A list of non-zero pixels can be returned using ``.nonzero_pixels()`` (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. + Map Interpolation ================= From d2cdd070cd91b84e0faa6e25be59d2935074873e Mon Sep 17 00:00:00 2001 From: Alexandra Rahlin Date: Tue, 2 Aug 2022 18:17:42 -0500 Subject: [PATCH 46/55] Mention multiplication, memory usage tricks --- maps/README.rst | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/maps/README.rst b/maps/README.rst index 155b68bd..a1838a57 100644 --- a/maps/README.rst +++ b/maps/README.rst @@ -67,7 +67,7 @@ Masks are returned when using comparison operators with map objects, e.g. ``map 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. A list of non-zero pixels can be returned using ``.nonzero_pixels()`` (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()``. +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_pixels()`` (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. @@ -76,6 +76,8 @@ 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. + Map Interpolation ================= From f0fd103dc3d501452cfffbfb3819b963c3dd498c Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Tue, 2 Aug 2022 17:34:51 -0500 Subject: [PATCH 47/55] tests cleanup --- maps/tests/flatsky_maps.py | 2 +- maps/tests/flatskymap_operators.py | 8 +++++++- maps/tests/healpix_maps.py | 2 +- 3 files changed, 9 insertions(+), 3 deletions(-) diff --git a/maps/tests/flatsky_maps.py b/maps/tests/flatsky_maps.py index 9c8eade9..ba0fcf19 100755 --- a/maps/tests/flatsky_maps.py +++ b/maps/tests/flatsky_maps.py @@ -122,4 +122,4 @@ assert not (set(pix1) ^ set(pix2)) mask = maps.make_point_source_mask(x, numpy.asarray(avec), numpy.asarray(dvec), numpy.asarray(rvec)) -assert not (set(numpy.where(numpy.asarray(mask.to_map()).ravel())[0]) ^ masked) +assert not (set(mask.nonzero_pixels()) ^ masked) diff --git a/maps/tests/flatskymap_operators.py b/maps/tests/flatskymap_operators.py index 63a5bb8d..576f9abc 100755 --- a/maps/tests/flatskymap_operators.py +++ b/maps/tests/flatskymap_operators.py @@ -264,12 +264,18 @@ med2 = get_map_median(mpad, ignore_zeros=True) assert(np.allclose(med2, med0)) - np.asarray(mpad)[np.asarray(mpad) == 0] = np.nan + 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) 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 = get_map_median(mpad, mask=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 3e926bfe..234737ce 100755 --- a/maps/tests/healpix_maps.py +++ b/maps/tests/healpix_maps.py @@ -98,7 +98,7 @@ assert not (set(pix1) ^ set(pix2)) mask = maps.make_point_source_mask(x, np.asarray(avec), np.asarray(dvec), np.asarray(rvec)) -assert not (set(np.where(np.asarray(mask.to_map()).ravel())[0]) ^ masked) +assert not (set(mask.nonzero_pixels()) ^ masked) x.shift_ra = True From 4f5c69ef43ba4f032c8e6248979a50f71ad39244 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Tue, 20 Sep 2022 13:36:53 -0500 Subject: [PATCH 48/55] create masks from numpy arrays via constructor or array_clone --- maps/include/maps/G3SkyMapMask.h | 8 ++ maps/src/G3SkyMapMask.cxx | 151 +++++++++++++++++++++++++++++-- maps/tests/mask_operators.py | 13 ++- 3 files changed, 161 insertions(+), 11 deletions(-) diff --git a/maps/include/maps/G3SkyMapMask.h b/maps/include/maps/G3SkyMapMask.h index cf222272..7d56fe43 100644 --- a/maps/include/maps/G3SkyMapMask.h +++ b/maps/include/maps/G3SkyMapMask.h @@ -23,11 +23,15 @@ 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); @@ -110,6 +114,10 @@ class G3SkyMapMask : public G3FrameObject { 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"); }; diff --git a/maps/src/G3SkyMapMask.cxx b/maps/src/G3SkyMapMask.cxx index e7d6eca5..194a103a 100644 --- a/maps/src/G3SkyMapMask.cxx +++ b/maps/src/G3SkyMapMask.cxx @@ -18,16 +18,31 @@ G3SkyMapMask::G3SkyMapMask(const G3SkyMap &parent, bool use_data, tmp->weighted = false; parent_ = tmp; data_ = std::vector(parent.size()); - if (use_data) { - for (size_t i = 0; i < parent.size(); i++) { - double v = parent.at(i); - if (zero_nans && v != v) - continue; - if (zero_infs && !std::isfinite(v)) - continue; - data_[i] = (v != 0); - } + + 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() @@ -45,6 +60,107 @@ G3SkyMapMask::Clone(bool copy_data) const 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 (zero_nans && v != v) + continue; + if (zero_infs && !std::isfinite(v)) + continue; + data_[i] = (v != 0); + } +} + +#define FILL_BUFFER(type) \ + for (size_t i = 0; i < view.len / sizeof(type); i++) { \ + double v = ((type *)view.buf)[i]; \ + if (zero_nans && v != v) \ + continue; \ + if (zero_infs && !std::isfinite(v)) \ + continue; \ + data_[i] = (v != 0); \ + } + +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) { @@ -404,15 +520,30 @@ PYBINDINGS("maps") { using namespace boost::python; - EXPORT_FRAMEOBJECT_NOINITNAMESPACE(G3SkyMapMask, (init((arg("parent"), arg("use_data")=false, arg("zero_nans")=false, arg("zero_infs")=false))), + 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.") diff --git a/maps/tests/mask_operators.py b/maps/tests/mask_operators.py index dd700158..9d65022e 100755 --- a/maps/tests/mask_operators.py +++ b/maps/tests/mask_operators.py @@ -5,7 +5,7 @@ maplist = [ maps.FlatSkyMap(np.random.randn(300, 500), core.G3Units.arcmin), - maps.HealpixSkyMap(np.random.randn(12 * 256 * 256)), + maps.HealpixSkyMap(np.random.randn(12 * 128 * 128)), ] for x in maplist: @@ -42,6 +42,17 @@ 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 m1 = x == 0 m2 = x != 0 From 962b4ee3e32645c67a8f368029734817be8a8195 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Tue, 20 Sep 2022 13:42:32 -0500 Subject: [PATCH 49/55] consistent logic for filling masks --- maps/src/G3SkyMapMask.cxx | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/maps/src/G3SkyMapMask.cxx b/maps/src/G3SkyMapMask.cxx index 194a103a..91b30c1a 100644 --- a/maps/src/G3SkyMapMask.cxx +++ b/maps/src/G3SkyMapMask.cxx @@ -67,22 +67,26 @@ G3SkyMapMask::FillFromMap(const G3SkyMap &m, bool zero_nans, bool zero_infs) for (size_t i = 0; i < m.size(); i++) { double v = m.at(i); - if (zero_nans && v != v) + if (v == 0) continue; - if (zero_infs && !std::isfinite(v)) + if (zero_nans && std::isnan(v)) continue; - data_[i] = (v != 0); + 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 (zero_nans && v != v) \ + if (v == 0) \ + continue; \ + if (zero_nans && std::isnan(v)) \ continue; \ - if (zero_infs && !std::isfinite(v)) \ + if (zero_infs && std::isinf(v)) \ continue; \ - data_[i] = (v != 0); \ + data_[i] = true; \ } void From 9ddc563c13be70b805c442d3382c008eeb7f6a6b Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Mon, 26 Sep 2022 17:41:20 -0500 Subject: [PATCH 50/55] Consistent interface for reduction ufuncs for maps and masks --- maps/include/maps/FlatSkyMap.h | 1 - maps/include/maps/G3SkyMap.h | 29 ++- maps/include/maps/HealpixSkyMap.h | 1 - maps/include/maps/maputils.h | 10 - maps/python/skymapaddons.py | 25 ++ maps/src/FlatSkyMap.cxx | 8 - maps/src/G3SkyMap.cxx | 377 ++++++++++++++++++++++++++++- maps/src/G3SkyMapMask.cxx | 19 +- maps/src/HealpixSkyMap.cxx | 8 - maps/src/maputils.cxx | 81 ------- maps/tests/flatsky_maps.py | 2 +- maps/tests/flatskymap_operators.py | 12 +- maps/tests/healpix_maps.py | 2 +- maps/tests/mask_operators.py | 67 ++++- 14 files changed, 505 insertions(+), 137 deletions(-) diff --git a/maps/include/maps/FlatSkyMap.h b/maps/include/maps/FlatSkyMap.h index 793c049f..7f97e50f 100644 --- a/maps/include/maps/FlatSkyMap.h +++ b/maps/include/maps/FlatSkyMap.h @@ -94,7 +94,6 @@ class FlatSkyMap : public G3FrameObject, public G3SkyMap { void NonZeroPixels(std::vector &indices, std::vector &data) const; void ApplyMask(const G3SkyMapMask &mask, bool inverse=false) override; - bool any() const override; void SetProj(MapProjection proj); void SetAlphaCenter(double alpha); diff --git a/maps/include/maps/G3SkyMap.h b/maps/include/maps/G3SkyMap.h index c29519e1..bf5e236d 100644 --- a/maps/include/maps/G3SkyMap.h +++ b/maps/include/maps/G3SkyMap.h @@ -124,8 +124,31 @@ class G3SkyMap { virtual G3SkyMapMask operator>(double rhs); virtual G3SkyMapMask operator>=(double rhs); - virtual bool all() const; - virtual bool any() const; + 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); @@ -176,7 +199,7 @@ class G3SkyMap { throw std::runtime_error("Compactification not implemented"); } - virtual boost::shared_ptr MakeMask(bool zero_nans = false, + virtual G3SkyMapMaskPtr MakeMask(bool zero_nans = false, bool zero_infs = false) const; protected: diff --git a/maps/include/maps/HealpixSkyMap.h b/maps/include/maps/HealpixSkyMap.h index fc75c359..b9ccfb31 100644 --- a/maps/include/maps/HealpixSkyMap.h +++ b/maps/include/maps/HealpixSkyMap.h @@ -72,7 +72,6 @@ class HealpixSkyMap : public G3FrameObject, public G3SkyMap { void NonZeroPixels(std::vector &indices, std::vector &data) const; // Iterators better? void ApplyMask(const G3SkyMapMask &mask, bool inverse=false) override; - bool any() const override; size_t nside() const {return info_.nside();} bool nested() const {return info_.nested();} diff --git a/maps/include/maps/maputils.h b/maps/include/maps/maputils.h index 2510de13..c8935ca8 100644 --- a/maps/include/maps/maputils.h +++ b/maps/include/maps/maputils.h @@ -39,16 +39,6 @@ void FlattenPol(FlatSkyMapPtr Q, FlatSkyMapPtr U, G3SkyMapWeightsPtr W=NULL, dou 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 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, G3SkyMapMaskConstPtr 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, G3SkyMapMaskConstPtr mask=NULL, - 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 diff --git a/maps/python/skymapaddons.py b/maps/python/skymapaddons.py index 125720c4..2f157232 100644 --- a/maps/python/skymapaddons.py +++ b/maps/python/skymapaddons.py @@ -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): diff --git a/maps/src/FlatSkyMap.cxx b/maps/src/FlatSkyMap.cxx index b1bce27c..2a6a1f99 100644 --- a/maps/src/FlatSkyMap.cxx +++ b/maps/src/FlatSkyMap.cxx @@ -664,14 +664,6 @@ size_t FlatSkyMap::NpixNonZero() const { return 0; } -bool FlatSkyMap::any() const -{ - for (auto i: *this) - if (i.second != 0) - return true; - return false; -} - #define GETSET(name, cname, type) \ type FlatSkyMap::name() const \ { \ diff --git a/maps/src/G3SkyMap.cxx b/maps/src/G3SkyMap.cxx index 70f0847d..2f6e58c6 100644 --- a/maps/src/G3SkyMap.cxx +++ b/maps/src/G3SkyMap.cxx @@ -6,6 +6,7 @@ #include #include #include +#include namespace bp=boost::python; @@ -667,23 +668,345 @@ skymap_compd(>=) skymap_compd(>) bool -G3SkyMap::all() const +G3SkyMap::all(G3SkyMapMaskConstPtr where) const { - for (size_t i = 0; i < size(); i++) - if (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; } bool -G3SkyMap::any() const +G3SkyMap::any(G3SkyMapMaskConstPtr where) const { - for (size_t i = 0; i < size(); i++) - if (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) { @@ -694,6 +1017,12 @@ 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) { @@ -1057,6 +1386,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")), @@ -1193,8 +1523,35 @@ PYBINDINGS("maps") { .def(bp::self > double()) .def("__bool__", &pyskymap_bool) - .def("any", &G3SkyMap::any) - .def("all", &G3SkyMap::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(); diff --git a/maps/src/G3SkyMapMask.cxx b/maps/src/G3SkyMapMask.cxx index 91b30c1a..ece15b9e 100644 --- a/maps/src/G3SkyMapMask.cxx +++ b/maps/src/G3SkyMapMask.cxx @@ -520,6 +520,16 @@ 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; @@ -556,11 +566,12 @@ PYBINDINGS("maps") .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("all", &G3SkyMapMask::all, "Test whether all elements are non-zero") - .def("any", &G3SkyMapMask::any, "Test whether any elements are non-zero") - .def("sum", &G3SkyMapMask::sum, "Sum of all elements in mask") - .def("nonzero_pixels", &G3SkyMapMask::NonZeroPixels, + .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), diff --git a/maps/src/HealpixSkyMap.cxx b/maps/src/HealpixSkyMap.cxx index f41c350d..cc52ea7d 100644 --- a/maps/src/HealpixSkyMap.cxx +++ b/maps/src/HealpixSkyMap.cxx @@ -910,14 +910,6 @@ size_t HealpixSkyMap::NpixNonZero() const return sz; } -bool HealpixSkyMap::any() const -{ - for (auto i: *this) - if (i.second != 0) - return true; - return false; -} - size_t HealpixSkyMap::AngleToPixel(double alpha, double delta) const { diff --git a/maps/src/maputils.cxx b/maps/src/maputils.cxx index 703a698f..4cd306bd 100644 --- a/maps/src/maputils.cxx +++ b/maps/src/maputils.cxx @@ -398,73 +398,6 @@ std::vector GetMapMoments(G3SkyMapConstPtr m, G3SkyMapMaskConstPtr mask, } -double -GetMapMedian(G3SkyMapConstPtr m, G3SkyMapMaskConstPtr mask, bool ignore_zeros, bool ignore_nans, - bool ignore_infs) -{ - - std::vector data; - size_t npix = mask ? mask->sum() : (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, G3SkyMapMaskConstPtr 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, G3SkyMapMaskConstPtr mask, bool ignore_zeros, bool ignore_nans, bool ignore_infs) @@ -657,20 +590,6 @@ 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")=G3SkyMapMaskConstPtr(), - 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")=G3SkyMapMaskConstPtr(), - 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")=G3SkyMapMaskConstPtr(), bp::arg("ignore_zeros")=false, bp::arg("ignore_nans")=false, bp::arg("ignore_infs")=false), diff --git a/maps/tests/flatsky_maps.py b/maps/tests/flatsky_maps.py index ba0fcf19..33d5358a 100755 --- a/maps/tests/flatsky_maps.py +++ b/maps/tests/flatsky_maps.py @@ -122,4 +122,4 @@ assert not (set(pix1) ^ set(pix2)) mask = maps.make_point_source_mask(x, numpy.asarray(avec), numpy.asarray(dvec), numpy.asarray(rvec)) -assert not (set(mask.nonzero_pixels()) ^ masked) +assert not (set(mask.nonzero()) ^ masked) diff --git a/maps/tests/flatskymap_operators.py b/maps/tests/flatskymap_operators.py index 576f9abc..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_moments, 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 @@ -255,25 +255,25 @@ stats0 = [np.mean(m1), np.var(m1), skew(m1), kurtosis(m1)] 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_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)) 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 = get_map_median(mpad, mask=mask) + med4 = mpad.median(where=mask) assert(np.allclose(med4, med0)) # convolution diff --git a/maps/tests/healpix_maps.py b/maps/tests/healpix_maps.py index 234737ce..50f8f48c 100755 --- a/maps/tests/healpix_maps.py +++ b/maps/tests/healpix_maps.py @@ -98,7 +98,7 @@ assert not (set(pix1) ^ set(pix2)) mask = maps.make_point_source_mask(x, np.asarray(avec), np.asarray(dvec), np.asarray(rvec)) -assert not (set(mask.nonzero_pixels()) ^ masked) +assert not (set(mask.nonzero()) ^ masked) x.shift_ra = True diff --git a/maps/tests/mask_operators.py b/maps/tests/mask_operators.py index 9d65022e..89ccd28c 100755 --- a/maps/tests/mask_operators.py +++ b/maps/tests/mask_operators.py @@ -4,11 +4,12 @@ from spt3g import core, maps maplist = [ - maps.FlatSkyMap(np.random.randn(300, 500), core.G3Units.arcmin), - maps.HealpixSkyMap(np.random.randn(12 * 128 * 128)), + 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 @@ -30,6 +31,7 @@ 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()) @@ -54,6 +56,7 @@ assert(((x > 0) == mp).all()) # float comparison + print("float comp") m1 = x == 0 m2 = x != 0 assert(m1.sum() + m2.sum() == x.size) @@ -67,6 +70,7 @@ assert(m1.sum() + m2.sum() == x.size) # map comparison + print("map comp") y = 2 * x.clone() m1 = y == x @@ -82,6 +86,7 @@ assert(m1.sum() + m2.sum() == x.size) # logic operators + print("logic") assert((m1 == ~m2).all()) assert(not (m1 & m2).any()) assert((m1 | m2).all()) @@ -101,6 +106,7 @@ assert((m3 == m2).all()) # convert to mask + print("conversion") x2 = x.clone() x2 *= m1 assert((x2.to_mask() == m1).all()) @@ -114,10 +120,65 @@ assert(np.sum(x4) == m1.sum()) assert((x4.to_mask() == m1).all()) - mpix = m1.nonzero_pixels() + 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 + v1 = getattr(np, attr)(x, **kwargs) + 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() From f5ea4c64f4fa4aa3589525c2bcfc41313f1daad0 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Mon, 26 Sep 2022 17:49:27 -0500 Subject: [PATCH 51/55] basically premature optimization --- maps/src/MapTODMasker.cxx | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/maps/src/MapTODMasker.cxx b/maps/src/MapTODMasker.cxx index d94036a3..0c956ca6 100644 --- a/maps/src/MapTODMasker.cxx +++ b/maps/src/MapTODMasker.cxx @@ -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,8 +107,8 @@ 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_->Parent()->coord_ref, alpha, delta); - auto detpointing = mask_->Parent()->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++) From 84b0b695bc5424decc19794d4d8d4857d1dc2798 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Mon, 26 Sep 2022 18:03:54 -0500 Subject: [PATCH 52/55] update docs --- maps/README.rst | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/maps/README.rst b/maps/README.rst index a1838a57..4285005a 100644 --- a/maps/README.rst +++ b/maps/README.rst @@ -56,7 +56,7 @@ 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 ======= @@ -67,7 +67,7 @@ Masks are returned when using comparison operators with map objects, e.g. ``map 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_pixels()`` (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()``. +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. @@ -78,6 +78,13 @@ The current implementation of masks is to use a dense ``std::vector`` as t 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, such as ``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 ================= From baf3711eabfc5fd79bc67a4af74b2cb976a6cc0d Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Mon, 26 Sep 2022 18:06:12 -0500 Subject: [PATCH 53/55] update docs --- maps/README.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/maps/README.rst b/maps/README.rst index 4285005a..bff5b23b 100644 --- a/maps/README.rst +++ b/maps/README.rst @@ -81,9 +81,9 @@ In general, when working with high-resolution maps of any sort, it is important Statistics ========== -Most ``numpy.ufunc``-like methods are defined for map objects, such as ``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. +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 values can be tested using ``.isnan(), .isinf(), .isfinite()`` methods as well; these return ``G3SkyMapMask`` objects. Map Interpolation ================= From 3b4ed29f0f56b53a982428807522dc1aee1c3559 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Tue, 27 Sep 2022 17:58:33 -0500 Subject: [PATCH 54/55] backward compatible test for older numpy versions --- maps/tests/mask_operators.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/maps/tests/mask_operators.py b/maps/tests/mask_operators.py index 89ccd28c..7ede196e 100755 --- a/maps/tests/mask_operators.py +++ b/maps/tests/mask_operators.py @@ -145,7 +145,13 @@ kwargs = {"where": w} if w is not None else {} if attr in ["std", "var"]: kwargs["ddof"] = 1 - v1 = getattr(np, attr)(x, **kwargs) + 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) From d1b1240f70e2048f152a9e313d0544dc696d7797 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Mon, 3 Oct 2022 15:19:35 -0500 Subject: [PATCH 55/55] ensure compatible mask --- maps/src/FlatSkyMap.cxx | 2 ++ maps/src/HealpixSkyMap.cxx | 2 ++ 2 files changed, 4 insertions(+) diff --git a/maps/src/FlatSkyMap.cxx b/maps/src/FlatSkyMap.cxx index 74fe2948..f56c98e1 100644 --- a/maps/src/FlatSkyMap.cxx +++ b/maps/src/FlatSkyMap.cxx @@ -642,6 +642,8 @@ 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; diff --git a/maps/src/HealpixSkyMap.cxx b/maps/src/HealpixSkyMap.cxx index 664ba75e..3d006770 100644 --- a/maps/src/HealpixSkyMap.cxx +++ b/maps/src/HealpixSkyMap.cxx @@ -947,6 +947,8 @@ 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;