From 7f2a23d5b1639e49547c27e4aebdc1a0d95cbd14 Mon Sep 17 00:00:00 2001 From: froejdh_e Date: Tue, 10 Dec 2024 22:00:12 +0100 Subject: [PATCH] accumulating clusters in one array --- include/aare/ClusterFinder.hpp | 44 ++++++++------ include/aare/ClusterVector.hpp | 107 ++++++++++++++++++++++++--------- python/examples/play.py | 37 +++++++----- python/src/cluster.hpp | 19 ++++-- 4 files changed, 137 insertions(+), 70 deletions(-) diff --git a/include/aare/ClusterFinder.hpp b/include/aare/ClusterFinder.hpp index 33a00ea..5bd8866 100644 --- a/include/aare/ClusterFinder.hpp +++ b/include/aare/ClusterFinder.hpp @@ -23,7 +23,7 @@ enum class eventType { UNDEFINED_EVENT = -1 /** undefined */ }; -template +template 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 m_pedestal; + ClusterVector m_clusters; + public: ClusterFinder(Shape<2> image_size, Shape<2> cluster_size, @@ -42,9 +44,10 @@ 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]) { - fmt::print("TypeIndex: {}\n", sizeof(Dtype)); - }; + m_pedestal(image_size[0], image_size[1]), + m_clusters(m_cluster_sizeX, m_cluster_sizeY, 1'000'000) { + // clusters = ClusterVector(m_cluster_sizeX, m_cluster_sizeY, 2000); + }; void push_pedestal_frame(NDView frame) { m_pedestal.push(frame); @@ -54,17 +57,20 @@ class ClusterFinder { NDArray noise() { return m_pedestal.std(); } - ClusterVector - find_clusters_without_threshold(NDView frame) { - // std::vector clusters; - // std::vector clusters; //Hard coded 3x3 cluster - // clusters.reserve(2000); - ClusterVector clusters(m_cluster_sizeX, m_cluster_sizeY); + ClusterVector steal_clusters() { + ClusterVector tmp = std::move(m_clusters); + m_clusters = ClusterVector(m_cluster_sizeX, m_cluster_sizeY, 2000); + return tmp; + } + void + find_clusters(NDView frame) { + // // size_t capacity = 2000; + // // ClusterVector clusters(m_cluster_sizeX, m_cluster_sizeY, capacity); eventType event_type = eventType::PEDESTAL; - // TODO! deal with even size clusters - // currently 3,3 -> +/- 1 - // 4,4 -> +/- 2 + // // TODO! deal with even size clusters + // // currently 3,3 -> +/- 1 + // // 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 cluster_data(m_cluster_sizeX * + std::vector 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 = - static_cast( + CT tmp = + static_cast( 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(cluster_data.data())); } } } - return clusters; + // return clusters; } // template diff --git a/include/aare/ClusterVector.hpp b/include/aare/ClusterVector.hpp index 6016e87..5445cbf 100644 --- a/include/aare/ClusterVector.hpp +++ b/include/aare/ClusterVector.hpp @@ -1,33 +1,74 @@ #pragma once #include #include +#include + +#include template 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) - : m_cluster_size_x(cluster_size_x), m_cluster_size_y(cluster_size_y) { - size_t num_bytes = element_offset() * m_capacity; - m_data = new std::byte[num_bytes]{}; - // fmt::print("Allocating {} bytes\n", num_bytes); + ClusterVector(int32_t cluster_size_x, int32_t cluster_size_y, + size_t capacity = 1024) + : m_cluster_size_x(cluster_size_x), m_cluster_size_y(cluster_size_y), + m_capacity(capacity) { + 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; - 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; + allocate_buffer(m_capacity * 2); } std::byte *ptr = element_ptr(m_size); *reinterpret_cast(ptr) = x; @@ -40,16 +81,17 @@ template class ClusterVector { m_size++; } - std::vector sum(){ + std::vector sum() { std::vector 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; - std::byte *ptr = element_ptr(i) + 2 * sizeof(int32_t); - for (size_t j = 0; j < m_cluster_size_x * m_cluster_size_y; j++) { - sum += *reinterpret_cast(ptr); - ptr += sizeof(T); - } - sums[i] = sum; + sums[i] = + std::accumulate(reinterpret_cast(ptr), + reinterpret_cast(ptr) + n_pixels, T{}); + ptr += stride; } return sums; } @@ -59,18 +101,25 @@ template 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 { - return element_offset() * i; - } + size_t element_offset(size_t i) const { return element_offset() * i; } - std::byte* element_ptr(size_t i) { - return m_data + element_offset(i); - } + std::byte *element_ptr(size_t 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; + } }; \ No newline at end of file diff --git a/python/examples/play.py b/python/examples/play.py index c5abcbf..986b718 100644 --- a/python/examples/play.py +++ b/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'Cluster per frame {total_clusters/N:.3f}') \ No newline at end of file +# print(f'Clustering {N} frames took {t_elapsed:.2f}s {N/t_elapsed:.0f} FPS') + + +# 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}') \ No newline at end of file diff --git a/python/src/cluster.hpp b/python/src/cluster.hpp index b9296fd..480aea1 100644 --- a/python/src/cluster.hpp +++ b/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::sum) + .def("sum", [](ClusterVector &self) { + auto *vec = new std::vector(self.sum()); + return return_vector(vec); + }) .def_buffer([typestr](ClusterVector &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 &self) { + auto v = new ClusterVector(self.steal_clusters()); + return v; + }) + .def("find_clusters", [](ClusterFinder &self, py::array_t frame) { auto view = make_view_2d(frame); - auto *vec = new ClusterVector( - self.find_clusters_without_threshold(view)); - return vec; + self.find_clusters(view); + return; }); m.def("hello", []() {