accumulating clusters in one array

include/aare/ClusterFinder.hpp

@@ -23,7 +23,7 @@ enum class eventType {
     UNDEFINED_EVENT = -1 /** undefined */
 };
 
-template <typename FRAME_TYPE = uint16_t, typename PEDESTAL_TYPE = double>
+template <typename FRAME_TYPE = uint16_t, typename PEDESTAL_TYPE = double, typename CT = int32_t>
 class ClusterFinder {
     Shape<2> m_image_size;
     const int m_cluster_sizeX;
@@ -33,6 +33,8 @@ class ClusterFinder {
     const double c2;
     const double c3;
     Pedestal<PEDESTAL_TYPE> m_pedestal;
+    ClusterVector<CT> m_clusters;
+
 
   public:
     ClusterFinder(Shape<2> image_size, Shape<2> cluster_size,
@@ -42,8 +44,9 @@ class ClusterFinder {
           m_nSigma(nSigma),
           c2(sqrt((m_cluster_sizeY + 1) / 2 * (m_cluster_sizeX + 1) / 2)),
           c3(sqrt(m_cluster_sizeX * m_cluster_sizeY)),
-          m_pedestal(image_size[0], image_size[1]) {
+          m_pedestal(image_size[0], image_size[1]),
-        fmt::print("TypeIndex: {}\n", sizeof(Dtype));
+          m_clusters(m_cluster_sizeX, m_cluster_sizeY, 1'000'000) {
+            // clusters = ClusterVector<CT>(m_cluster_sizeX, m_cluster_sizeY, 2000);
           };
 
     void push_pedestal_frame(NDView<FRAME_TYPE, 2> frame) {
@@ -54,17 +57,20 @@ class ClusterFinder {
 
     NDArray<PEDESTAL_TYPE, 2> noise() { return m_pedestal.std(); }
 
-    ClusterVector<PEDESTAL_TYPE>
+    ClusterVector<CT> steal_clusters() {
-    find_clusters_without_threshold(NDView<FRAME_TYPE, 2> frame) {
+        ClusterVector<CT> tmp = std::move(m_clusters);
-        // std::vector<DynamicCluster> clusters;
+        m_clusters = ClusterVector<CT>(m_cluster_sizeX, m_cluster_sizeY, 2000);
-        // std::vector<Cluster> clusters; //Hard coded 3x3 cluster
+        return tmp;
-        // clusters.reserve(2000);
+    }
-        ClusterVector<PEDESTAL_TYPE> clusters(m_cluster_sizeX, m_cluster_sizeY);
+    void
+    find_clusters(NDView<FRAME_TYPE, 2> frame) {
+        // // size_t capacity = 2000;
+        // // ClusterVector<CT> clusters(m_cluster_sizeX, m_cluster_sizeY, capacity);
         eventType event_type = eventType::PEDESTAL;
 
-        // TODO! deal with even size clusters
+        // // TODO! deal with even size clusters
-        // currently 3,3 -> +/- 1
+        // // currently 3,3 -> +/- 1
-        //  4,4 -> +/- 2
+        // //  4,4 -> +/- 2
         short dy = m_cluster_sizeY / 2;
         short dx = m_cluster_sizeX / 2;
 
@@ -108,29 +114,29 @@ class ClusterFinder {
                     event_type = eventType::PHOTON_MAX;
 
                     short i = 0;
-                    std::vector<PEDESTAL_TYPE> cluster_data(m_cluster_sizeX *
+                    std::vector<CT> cluster_data(m_cluster_sizeX *
                                                             m_cluster_sizeY);
 
                     for (short ir = -dy; ir < dy + 1; ir++) {
                         for (short ic = -dx; ic < dx + 1; ic++) {
                             if (ix + ic >= 0 && ix + ic < frame.shape(1) &&
                                 iy + ir >= 0 && iy + ir < frame.shape(0)) {
-                                PEDESTAL_TYPE tmp =
+                                CT tmp =
-                                    static_cast<PEDESTAL_TYPE>(
+                                    static_cast<CT>(
                                         frame(iy + ir, ix + ic)) -
                                     m_pedestal.mean(iy + ir, ix + ic);
-                                cluster_data[i] = tmp;
+                                cluster_data[i] = tmp; //Watch for out of bounds access
                                 i++;
                             }
                         }
                     }
-                    clusters.push_back(
+                    m_clusters.push_back(
                         ix, iy,
                         reinterpret_cast<std::byte *>(cluster_data.data()));
                 }
             }
         }
-        return clusters;
+        // return clusters;
     }
 
     // template <typename FRAME_TYPE, typename PEDESTAL_TYPE>

include/aare/ClusterVector.hpp

@@ -1,33 +1,74 @@
 #pragma once
 #include <cstddef>
 #include <cstdint>
+#include <numeric>
+
+#include <fmt/core.h>
 
 template <typename T> class ClusterVector {
     int32_t m_cluster_size_x;
     int32_t m_cluster_size_y;
     std::byte *m_data{};
     size_t m_size{0};
-    size_t m_capacity{10};
+    size_t m_capacity;
 
   public:
-    ClusterVector(int32_t cluster_size_x, int32_t cluster_size_y)
+    ClusterVector(int32_t cluster_size_x, int32_t cluster_size_y,
-        : m_cluster_size_x(cluster_size_x), m_cluster_size_y(cluster_size_y) {
+                  size_t capacity = 1024)
-        size_t num_bytes = element_offset() * m_capacity;
+        : m_cluster_size_x(cluster_size_x), m_cluster_size_y(cluster_size_y),
-        m_data = new std::byte[num_bytes]{};
+          m_capacity(capacity) {
-        // fmt::print("Allocating {} bytes\n", num_bytes);
+        allocate_buffer(capacity);
+    }
+    ~ClusterVector() {
+        fmt::print("~ClusterVector - size: {}, capacity: {}\n", m_size,
+                   m_capacity);
+        delete[] m_data;
+    }
+
+    // ClusterVector(const ClusterVector & other){
+    //     m_cluster_size_x = other.m_cluster_size_x;
+    //     m_cluster_size_y = other.m_cluster_size_y;
+    //     m_data = new std::byte[other.m_capacity];
+    //     std::copy(other.m_data, other.m_data + other.m_capacity, m_data);
+    //     m_size = other.m_size;
+    //     m_capacity = other.m_capacity;
+    // }
+
+    ClusterVector(ClusterVector &&other) noexcept
+        : m_cluster_size_x(other.m_cluster_size_x),
+          m_cluster_size_y(other.m_cluster_size_y), m_data(other.m_data),
+          m_size(other.m_size), m_capacity(other.m_capacity) {
+        other.m_data = nullptr;
+        other.m_size = 0;
+        other.m_capacity = 0;
+    }
+
+    //Move assignment operator
+    ClusterVector& operator=(ClusterVector &&other) noexcept {
+        if (this != &other) {
+            delete[] m_data;
+            m_cluster_size_x = other.m_cluster_size_x;
+            m_cluster_size_y = other.m_cluster_size_y;
+            m_data = other.m_data;
+            m_size = other.m_size;
+            m_capacity = other.m_capacity;
+            other.m_data = nullptr;
+            other.m_size = 0;
+            other.m_capacity = 0;
+        }
+        return *this;
+    }
+
+    void reserve(size_t capacity) {
+        if (capacity > m_capacity) {
+            allocate_buffer(capacity);
+        }
     }
 
     // data better hold data of the right size!
     void push_back(int32_t x, int32_t y, const std::byte *data) {
         if (m_size == m_capacity) {
-            m_capacity *= 2;
+            allocate_buffer(m_capacity * 2);
-            std::byte *new_data =
-                new std::byte[element_offset()*m_capacity]{};
-            std::copy(m_data,
-                      m_data + element_offset()*m_size,
-                      new_data);
-            delete[] m_data;
-            m_data = new_data;
         }
         std::byte *ptr = element_ptr(m_size);
         *reinterpret_cast<int32_t *>(ptr) = x;
@@ -42,14 +83,15 @@ template <typename T> class ClusterVector {
 
     std::vector<T> sum() {
         std::vector<T> sums(m_size);
+        const size_t stride = element_offset();
+        const size_t n_pixels = m_cluster_size_x * m_cluster_size_y;
+        std::byte *ptr = m_data + 2 * sizeof(int32_t); // skip x and y
+
         for (size_t i = 0; i < m_size; i++) {
-            T sum = 0;
+            sums[i] =
-            std::byte *ptr = element_ptr(i) + 2 * sizeof(int32_t);
+                std::accumulate(reinterpret_cast<T *>(ptr),
-            for (size_t j = 0; j < m_cluster_size_x * m_cluster_size_y; j++) {
+                                reinterpret_cast<T *>(ptr) + n_pixels, T{});
-                sum += *reinterpret_cast<T *>(ptr);
+            ptr += stride;
-                ptr += sizeof(T);
-            }
-            sums[i] = sum;
         }
         return sums;
     }
@@ -59,18 +101,25 @@ template <typename T> class ClusterVector {
         return sizeof(m_cluster_size_x) + sizeof(m_cluster_size_y) +
                m_cluster_size_x * m_cluster_size_y * sizeof(T);
     }
-    size_t element_offset(size_t i) const {
+    size_t element_offset(size_t i) const { return element_offset() * i; }
-        return element_offset() * i;
-    }
 
-    std::byte* element_ptr(size_t i) {
+    std::byte *element_ptr(size_t i) { return m_data + element_offset(i); }
-        return m_data + element_offset(i);
-    }
 
     int16_t cluster_size_x() const { return m_cluster_size_x; }
     int16_t cluster_size_y() const { return m_cluster_size_y; }
 
     std::byte *data() { return m_data; }
 
-    ~ClusterVector() { delete[] m_data; }
+  private:
+    void allocate_buffer(size_t new_capacity) {
+        size_t num_bytes = element_offset() * new_capacity;
+        fmt::print(
+            "ClusterVector allocating {} elements for a total of {} bytes\n",
+            new_capacity, num_bytes);
+        std::byte *new_data = new std::byte[num_bytes]{};
+        std::copy(m_data, m_data + element_offset() * m_size, new_data);
+        delete[] m_data;
+        m_data = new_data;
+        m_capacity = new_capacity;
+    }
 };

python/examples/play.py

@@ -22,7 +22,9 @@ im = ax.imshow(cf.pedestal())
 cf.pedestal()
 cf.noise()
 
-N = 200
+
+
+N = 500
 t0 = time.perf_counter()
 hist1 = bh.Histogram(bh.axis.Regular(40, -2, 4000))
 f.seek(0)
@@ -31,23 +33,26 @@ t0 = time.perf_counter()
 data = f.read_n(N)
 t_elapsed = time.perf_counter()-t0
 
-print(f'Reading {N} frames took {t_elapsed:.3f}s {N/t_elapsed:.0f} FPS')
 
-clusters = []
+n_bytes = data.itemsize*data.size
+
+print(f'Reading {N} frames took {t_elapsed:.3f}s {N/t_elapsed:.0f} FPS, {n_bytes/1024**2:.4f} GB/s')
+
+
 for frame in data:
-    clusters += [cf.find_clusters_without_threshold(frame)]
+    a = cf.find_clusters(frame)
 
+clusters = cf.steal_clusters()
 
-t_elapsed = time.perf_counter()-t0
-print(f'Clustering {N} frames took {t_elapsed:.2f}s  {N/t_elapsed:.0f} FPS')
-
-
-t0 = time.perf_counter()
-total_clusters = 0
-for cl in clusters:
-        arr = np.array(cl, copy = False)
-        hist1.fill(arr['data'].sum(axis = 1).sum(axis = 1))
-        total_clusters += cl.size
 # t_elapsed = time.perf_counter()-t0
-print(f'Filling histogram with {total_clusters} clusters took: {t_elapsed:.3f}s')
+# print(f'Clustering {N} frames took {t_elapsed:.2f}s  {N/t_elapsed:.0f} FPS')
-print(f'Cluster per frame {total_clusters/N:.3f}')
+
+
+# t0 = time.perf_counter()
+# total_clusters = clusters.size
+
+# hist1.fill(clusters.sum())
+
+# t_elapsed = time.perf_counter()-t0
+# print(f'Filling histogram with the sum of {total_clusters} clusters took: {t_elapsed:.3f}s, {total_clusters/t_elapsed:.3g} clust/s')
+# print(f'Average number of clusters per frame {total_clusters/N:.3f}')

python/src/cluster.hpp

@@ -26,12 +26,15 @@ void define_cluster_vector(py::module &m, const std::string &typestr) {
                      "i:x:\ni:y:\n({},{}){}:data:", v.cluster_size_x(),
                      v.cluster_size_y(), typestr);
              })
-        .def("sum", &ClusterVector<T>::sum)
+        .def("sum", [](ClusterVector<T> &self) {
+            auto *vec = new std::vector<T>(self.sum());
+            return return_vector(vec);
+        })
         .def_buffer([typestr](ClusterVector<T> &v) -> py::buffer_info {
             return py::buffer_info(
                 v.data(),           /* Pointer to buffer */
                 v.element_offset(), /* Size of one scalar */
-                fmt::format("i:x:\ni:y:\n({},{}){}:data:", v.cluster_size_x(),
+                fmt::format("i:x:\ni:y:\n{}{}:data:", v.cluster_size_x()*
                             v.cluster_size_y(),
                             typestr), /* Format descriptor */
                 1,                    /* Number of dimensions */
@@ -62,13 +65,17 @@ void define_cluster_finder_bindings(py::module &m) {
                  *arr = self.noise();
                  return return_image_data(arr);
              })
-        .def("find_clusters_without_threshold",
+        .def("steal_clusters",
+             [](ClusterFinder<uint16_t, pd_type> &self) {
+                 auto v = new ClusterVector<int>(self.steal_clusters());
+                 return v;
+             })
+        .def("find_clusters",
              [](ClusterFinder<uint16_t, pd_type> &self,
                 py::array_t<uint16_t> frame) {
                  auto view = make_view_2d(frame);
-                 auto *vec = new ClusterVector<pd_type>(
+                 self.find_clusters(view);
-                     self.find_clusters_without_threshold(view));
+                 return;
-                 return vec;
              });
 
     m.def("hello", []() {