Port MapSampler class and utilities.

src/libtoast/include/toast.hpp

@@ -16,7 +16,7 @@
 #include <toast/math_healpix.hpp>
 #include <toast/fod_psd.hpp>
 #include <toast/map_cov.hpp>
-#include <toast/todmap_scanning.hpp>
+#include <toast/tod_mapscan.hpp>
 #include <toast/tod_filter.hpp>
 #include <toast/tod_pointing.hpp>
 #include <toast/tod_simnoise.hpp>

src/libtoast/include/toast/todmap_scanning.hpp → src/libtoast/include/toast/tod_mapscan.hpp

@@ -3,11 +3,12 @@
 // All rights reserved.  Use of this source code is governed by
 // a BSD-style license that can be found in the LICENSE file.
 
-#ifndef TOAST_TODMAP_SCANNING_HPP
-#define TOAST_TODMAP_SCANNING_HPP
+#ifndef TOAST_TOD_MAPSCAN_HPP
+#define TOAST_TOD_MAPSCAN_HPP
 
-namespace toast {
+#include <cstring>
 
+namespace toast {
 template <typename T>
 void scan_local_map(int64_t const * submap, int64_t subnpix, double const * weights,
                     int64_t nmap, int64_t * subpix, T const * map, double * tod,
@@ -33,6 +34,25 @@ void scan_local_map(int64_t const * submap, int64_t subnpix, double const * weig
     return;
 }
 
+template <typename T>
+void fast_scanning(double * toi, int64_t nsamp,
+                   int64_t const * pixels, double const * weights,
+                   int64_t nweight, T const * bmap) {
+    memset(toi, 0, nsamp * sizeof(double));
+    #pragma omp parallel for
+    for (int64_t row = 0; row < nsamp; ++row) {
+        int64_t offset = row * nweight;
+        for (int64_t col = 0; col < nweight; ++col) {
+            int64_t pix = pixels[offset];
+            if (pix < 0) continue;
+            double weight = weights[offset];
+            toi[row] += bmap[pix] * weight;
+            ++offset;
+        }
+    }
+    return;
+}
+
 //
 // template <typename T>
 // void scan_global_map(int64_t npixmap, int64_t * pixels,
@@ -101,28 +121,6 @@ void scan_local_map(int64_t const * submap, int64_t subnpix, double const * weig
 //     return;
 // }
 //
-// void toast::map_tools::fast_scanning32(double * toi, int64_t const nsamp,
-//                                        int64_t const * pixels,
-//                                        double const * weights,
-//                                        int64_t const nweight,
-//                                        float const * bmap) {
-//     memset(toi, 0, nsamp * sizeof(double));
-//     #pragma \
-//     omp parallel for schedule(static) default(none) shared(pixels, weights, nweight,
-// map, tod, nsamp
-//     for (int64_t row = 0; row < nsamp; ++row) {
-//         int64_t offset = row * nweight;
-//         for (int64_t col = 0; col < nweight; ++col) {
-//             int64_t pix = pixels[offset];
-//             if (pix < 0) continue;
-//             double weight = weights[offset];
-//             toi[row] += bmap[pix] * weight;
-//             ++offset;
-//         }
-//     }
-// }
-//
-
 }
 
-#endif // ifndef TOAST_TODMAP_SCANNING_HPP
+#endif // ifndef TOAST_TOD_MAPSCAN_HPP

src/toast/_libtoast_todmap_scanning.cpp

@@ -74,10 +74,59 @@ void register_scan_map(py::module & m, char const * name) {
     return;
 }
 
+template <typename T>
+void register_fast_scanning(py::module & m, char const * name) {
+    m.def(name,
+          [](py::buffer tod, py::buffer pix, py::buffer weights, py::buffer mapdata) {
+              pybuffer_check_1D <int64_t> (pix);
+              pybuffer_check_1D <T> (mapdata);
+              pybuffer_check_1D <double> (weights);
+              pybuffer_check_1D <double> (tod);
+              py::buffer_info info_pix = pix.request();
+              py::buffer_info info_mapdata = mapdata.request();
+              py::buffer_info info_weights = weights.request();
+              py::buffer_info info_tod = tod.request();
+              size_t nsamp = info_tod.size;
+              if (info_pix.size != nsamp) {
+                  auto log = toast::Logger::get();
+                  std::ostringstream o;
+                  o << "Buffer sizes are not consistent.";
+                  log.error(o.str().c_str());
+                  throw std::runtime_error(o.str().c_str());
+              }
+              size_t nw = (size_t)(info_weights.size / nsamp);
+              int64_t * rawpix = reinterpret_cast <int64_t *> (info_pix.ptr);
+              T * rawmapdata = reinterpret_cast <T *> (info_mapdata.ptr);
+              double * rawweights = reinterpret_cast <double *> (info_weights.ptr);
+              double * rawtod = reinterpret_cast <double *> (info_tod.ptr);
+              toast::fast_scanning <T> (rawtod, nsamp, rawpix, rawweights,
+                                        nw, rawmapdata);
+              return;
+          }, py::arg("tod"), py::arg("pix"), py::arg("weights"), py::arg(
+              "mapdata"), R"(
+        Scan global maps into timestreams.
+
+        Args:
+            tod (array, float64):  The timestream on which to accumulate the map
+                values.
+            pix (array, int64):  For each time domain sample, the pixel index.
+            weights (array, float64):  The pointing matrix weights for each time
+                sample and map.
+            mapdata (array):  The flattened local piece of the map.
+
+        Returns:
+            None.
+
+    )");
+    return;
+}
+
 void init_todmap_scanning(py::module & m) {
     register_scan_map <double> (m, "scan_map_float64");
     register_scan_map <float> (m, "scan_map_float32");
     register_scan_map <int64_t> (m, "scan_map_int64");
     register_scan_map <int32_t> (m, "scan_map_int32");
+    register_fast_scanning <double> (m, "fast_scanning_float64");
+    register_fast_scanning <float> (m, "fast_scanning_float32");
     return;
 }

src/toast/map/CMakeLists.txt

@@ -4,7 +4,6 @@
 install(FILES
     __init__.py
     madam.py
-    mapsampler.py
     cov.py
     pixels.py
     pysm.py

src/toast/map/__init__.py

@@ -7,6 +7,7 @@
 from .pixels import OpLocalPixels, DistPixels
 
 from .pysm import pysm
+
 if pysm is not None:
     from .pysm import PySMSky
 
@@ -18,6 +19,4 @@
     covariance_apply,
 )
 
-# from .mapsampler import MapSampler
-
 from .madam import OpMadam