From b3a9e9576b5733e92498dc7f7b036d4e646c6e3b Mon Sep 17 00:00:00 2001 From: froejdh_e Date: Wed, 11 Dec 2024 16:27:36 +0100 Subject: [PATCH] WIP --- CMakeLists.txt | 1 + docs/conf.py.in | 1 - docs/src/ClusterVector.rst | 6 + docs/src/index.rst | 1 + include/aare/ClusterFinder.hpp | 305 +++++++++++++++++---------------- include/aare/ClusterVector.hpp | 52 +++++- python/src/cluster.hpp | 41 +---- src/ClusterVector.test.cpp | 77 +++++++++ 8 files changed, 301 insertions(+), 183 deletions(-) create mode 100644 docs/src/ClusterVector.rst create mode 100644 src/ClusterVector.test.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index f67d655..a38e8fd 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -344,6 +344,7 @@ if(AARE_TESTS) ${CMAKE_CURRENT_SOURCE_DIR}/src/NDArray.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/NDView.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/ClusterFinder.test.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/src/ClusterVector.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/Pedestal.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyFile.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyHelpers.test.cpp diff --git a/docs/conf.py.in b/docs/conf.py.in index 3702330..ad73575 100644 --- a/docs/conf.py.in +++ b/docs/conf.py.in @@ -29,7 +29,6 @@ version = '@PROJECT_VERSION@' # extensions coming with Sphinx (named 'sphinx.ext.*') or your custom # ones. extensions = ['breathe', - 'sphinx_rtd_theme', 'sphinx.ext.autodoc', 'sphinx.ext.napoleon', ] diff --git a/docs/src/ClusterVector.rst b/docs/src/ClusterVector.rst new file mode 100644 index 0000000..bb2a0d8 --- /dev/null +++ b/docs/src/ClusterVector.rst @@ -0,0 +1,6 @@ +ClusterVector +============= + +.. doxygenclass:: aare::ClusterVector + :members: + :undoc-members: \ No newline at end of file diff --git a/docs/src/index.rst b/docs/src/index.rst index 228d7c4..4316a2c 100644 --- a/docs/src/index.rst +++ b/docs/src/index.rst @@ -46,6 +46,7 @@ AARE Dtype ClusterFinder ClusterFile + ClusterVector Pedestal RawFile RawSubFile diff --git a/include/aare/ClusterFinder.hpp b/include/aare/ClusterFinder.hpp index 5bd8866..7111cf9 100644 --- a/include/aare/ClusterFinder.hpp +++ b/include/aare/ClusterFinder.hpp @@ -23,62 +23,83 @@ enum class eventType { UNDEFINED_EVENT = -1 /** undefined */ }; -template +template class ClusterFinder { Shape<2> m_image_size; const int m_cluster_sizeX; const int m_cluster_sizeY; - const double m_threshold; - const double m_nSigma; - const double c2; - const double c3; + // const PEDESTAL_TYPE m_threshold; + const PEDESTAL_TYPE m_nSigma; + const PEDESTAL_TYPE c2; + const PEDESTAL_TYPE c3; Pedestal m_pedestal; ClusterVector m_clusters; - public: + /** + * @brief Construct a new ClusterFinder object + * @param image_size size of the image + * @param cluster_size size of the cluster (x, y) + * @param nSigma number of sigma above the pedestal to consider a photon + * @param capacity initial capacity of the cluster vector + * + */ ClusterFinder(Shape<2> image_size, Shape<2> cluster_size, - double nSigma = 5.0, double threshold = 0.0) + PEDESTAL_TYPE nSigma = 5.0, size_t capacity = 1000000) : m_image_size(image_size), m_cluster_sizeX(cluster_size[0]), - m_cluster_sizeY(cluster_size[1]), m_threshold(threshold), + m_cluster_sizeY(cluster_size[1]), 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_clusters(m_cluster_sizeX, m_cluster_sizeY, 1'000'000) { - // clusters = ClusterVector(m_cluster_sizeX, m_cluster_sizeY, 2000); - }; + m_clusters(m_cluster_sizeX, m_cluster_sizeY, capacity) {}; void push_pedestal_frame(NDView frame) { m_pedestal.push(frame); } NDArray pedestal() { return m_pedestal.mean(); } - NDArray noise() { return m_pedestal.std(); } - ClusterVector steal_clusters() { + /** + * @brief Move the clusters from the ClusterVector in the ClusterFinder to a + * new ClusterVector and return it. + * @param realloc_same_capacity if true the new ClusterVector will have the + * same capacity as the old one + * + */ + ClusterVector steal_clusters(bool realloc_same_capacity = false) { ClusterVector tmp = std::move(m_clusters); - m_clusters = ClusterVector(m_cluster_sizeX, m_cluster_sizeY, 2000); + if (realloc_same_capacity) + m_clusters = ClusterVector(m_cluster_sizeX, m_cluster_sizeY, + tmp.capacity()); + else + m_clusters = ClusterVector(m_cluster_sizeX, m_cluster_sizeY); 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; - + void find_clusters(NDView frame) { // // 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; + std::vector cluster_data(m_cluster_sizeX * m_cluster_sizeY); for (int iy = 0; iy < frame.shape(0); iy++) { for (int ix = 0; ix < frame.shape(1); ix++) { + PEDESTAL_TYPE max = std::numeric_limits::min(); PEDESTAL_TYPE total = 0; + // What can we short circuit here? + PEDESTAL_TYPE rms = m_pedestal.std(iy, ix); + PEDESTAL_TYPE value = (frame(iy, ix) - m_pedestal.mean(iy, ix)); + + if (value < -m_nSigma * rms) + continue; // NEGATIVE_PEDESTAL go to next pixel + // TODO! No pedestal update??? + 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) && @@ -92,159 +113,157 @@ class ClusterFinder { } } } - PEDESTAL_TYPE rms = m_pedestal.std(iy, ix); - PEDESTAL_TYPE value = (frame(iy, ix) - m_pedestal.mean(iy, ix)); - if (value < -m_nSigma * rms) { - continue; // NEGATIVE_PEDESTAL go to next pixel - // TODO! No pedestal update??? - } else if (max > m_nSigma * rms) { - event_type = eventType::PHOTON; + if ((max > m_nSigma * rms)) { if (value < max) continue; // Not max go to the next pixel + // but also no pedestal update } else if (total > c3 * m_nSigma * rms) { - event_type = eventType::PHOTON; + // pass } else { m_pedestal.push(iy, ix, frame(iy, ix)); continue; // It was a pedestal value nothing to store } // Store cluster - if (event_type == eventType::PHOTON && value >= max) { - event_type = eventType::PHOTON_MAX; + if (value == max) { + // Zero out the cluster data + std::fill(cluster_data.begin(), cluster_data.end(), 0); - short i = 0; - 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++) { + // Fill the cluster data since we have a photon to store + // It's worth redoing the look since most of the time we + // don't have a photon + int i = 0; + for (int ir = -dy; ir < dy + 1; ir++) { + for (int ic = -dx; ic < dx + 1; ic++) { if (ix + ic >= 0 && ix + ic < frame.shape(1) && iy + ir >= 0 && iy + ir < frame.shape(0)) { CT tmp = - static_cast( - frame(iy + ir, ix + ic)) - + static_cast(frame(iy + ir, ix + ic)) - m_pedestal.mean(iy + ir, ix + ic); - cluster_data[i] = tmp; //Watch for out of bounds access + cluster_data[i] = + tmp; // Watch for out of bounds access i++; } } } + + // Add the cluster to the output ClusterVector m_clusters.push_back( ix, iy, reinterpret_cast(cluster_data.data())); } } } - // return clusters; } - // template - std::vector - find_clusters_with_threshold(NDView frame, - Pedestal &pedestal) { - assert(m_threshold > 0); - std::vector clusters; - std::vector> eventMask; - for (int i = 0; i < frame.shape(0); i++) { - eventMask.push_back(std::vector(frame.shape(1))); - } - double tthr, tthr1, tthr2; + // // template + // std::vector + // find_clusters_with_threshold(NDView frame, + // Pedestal &pedestal) { + // assert(m_threshold > 0); + // std::vector clusters; + // std::vector> eventMask; + // for (int i = 0; i < frame.shape(0); i++) { + // eventMask.push_back(std::vector(frame.shape(1))); + // } + // double tthr, tthr1, tthr2; - NDArray rest({frame.shape(0), frame.shape(1)}); - NDArray nph({frame.shape(0), frame.shape(1)}); - // convert to n photons - // nph = (frame-pedestal.mean()+0.5*m_threshold)/m_threshold; // can be - // optimized with expression templates? - for (int iy = 0; iy < frame.shape(0); iy++) { - for (int ix = 0; ix < frame.shape(1); ix++) { - auto val = frame(iy, ix) - pedestal.mean(iy, ix); - nph(iy, ix) = (val + 0.5 * m_threshold) / m_threshold; - nph(iy, ix) = nph(iy, ix) < 0 ? 0 : nph(iy, ix); - rest(iy, ix) = val - nph(iy, ix) * m_threshold; - } - } - // iterate over frame pixels - for (int iy = 0; iy < frame.shape(0); iy++) { - for (int ix = 0; ix < frame.shape(1); ix++) { - eventMask[iy][ix] = eventType::PEDESTAL; - // initialize max and total - FRAME_TYPE max = std::numeric_limits::min(); - long double total = 0; - if (rest(iy, ix) <= 0.25 * m_threshold) { - pedestal.push(iy, ix, frame(iy, ix)); - continue; - } - eventMask[iy][ix] = eventType::NEIGHBOUR; - // iterate over cluster pixels around the current pixel (ix,iy) - for (short ir = -(m_cluster_sizeY / 2); - ir < (m_cluster_sizeY / 2) + 1; ir++) { - for (short ic = -(m_cluster_sizeX / 2); - ic < (m_cluster_sizeX / 2) + 1; ic++) { - if (ix + ic >= 0 && ix + ic < frame.shape(1) && - iy + ir >= 0 && iy + ir < frame.shape(0)) { - auto val = frame(iy + ir, ix + ic) - - pedestal.mean(iy + ir, ix + ic); - total += val; - if (val > max) { - max = val; - } - } - } - } + // NDArray rest({frame.shape(0), frame.shape(1)}); + // NDArray nph({frame.shape(0), frame.shape(1)}); + // // convert to n photons + // // nph = (frame-pedestal.mean()+0.5*m_threshold)/m_threshold; // can + // be + // // optimized with expression templates? + // for (int iy = 0; iy < frame.shape(0); iy++) { + // for (int ix = 0; ix < frame.shape(1); ix++) { + // auto val = frame(iy, ix) - pedestal.mean(iy, ix); + // nph(iy, ix) = (val + 0.5 * m_threshold) / m_threshold; + // nph(iy, ix) = nph(iy, ix) < 0 ? 0 : nph(iy, ix); + // rest(iy, ix) = val - nph(iy, ix) * m_threshold; + // } + // } + // // iterate over frame pixels + // for (int iy = 0; iy < frame.shape(0); iy++) { + // for (int ix = 0; ix < frame.shape(1); ix++) { + // eventMask[iy][ix] = eventType::PEDESTAL; + // // initialize max and total + // FRAME_TYPE max = std::numeric_limits::min(); + // long double total = 0; + // if (rest(iy, ix) <= 0.25 * m_threshold) { + // pedestal.push(iy, ix, frame(iy, ix)); + // continue; + // } + // eventMask[iy][ix] = eventType::NEIGHBOUR; + // // iterate over cluster pixels around the current pixel + // (ix,iy) for (short ir = -(m_cluster_sizeY / 2); + // ir < (m_cluster_sizeY / 2) + 1; ir++) { + // for (short ic = -(m_cluster_sizeX / 2); + // ic < (m_cluster_sizeX / 2) + 1; ic++) { + // if (ix + ic >= 0 && ix + ic < frame.shape(1) && + // iy + ir >= 0 && iy + ir < frame.shape(0)) { + // auto val = frame(iy + ir, ix + ic) - + // pedestal.mean(iy + ir, ix + ic); + // total += val; + // if (val > max) { + // max = val; + // } + // } + // } + // } - auto rms = pedestal.std(iy, ix); - if (m_nSigma == 0) { - tthr = m_threshold; - tthr1 = m_threshold; - tthr2 = m_threshold; - } else { - tthr = m_nSigma * rms; - tthr1 = m_nSigma * rms * c3; - tthr2 = m_nSigma * rms * c2; + // auto rms = pedestal.std(iy, ix); + // if (m_nSigma == 0) { + // tthr = m_threshold; + // tthr1 = m_threshold; + // tthr2 = m_threshold; + // } else { + // tthr = m_nSigma * rms; + // tthr1 = m_nSigma * rms * c3; + // tthr2 = m_nSigma * rms * c2; - if (m_threshold > 2 * tthr) - tthr = m_threshold - tthr; - if (m_threshold > 2 * tthr1) - tthr1 = tthr - tthr1; - if (m_threshold > 2 * tthr2) - tthr2 = tthr - tthr2; - } - if (total > tthr1 || max > tthr) { - eventMask[iy][ix] = eventType::PHOTON; - nph(iy, ix) += 1; - rest(iy, ix) -= m_threshold; - } else { - pedestal.push(iy, ix, frame(iy, ix)); - continue; - } - if (eventMask[iy][ix] == eventType::PHOTON && - frame(iy, ix) - pedestal.mean(iy, ix) >= max) { - eventMask[iy][ix] = eventType::PHOTON_MAX; - DynamicCluster cluster(m_cluster_sizeX, m_cluster_sizeY, - Dtype(typeid(FRAME_TYPE))); - cluster.x = ix; - cluster.y = iy; - short i = 0; - for (short ir = -(m_cluster_sizeY / 2); - ir < (m_cluster_sizeY / 2) + 1; ir++) { - for (short ic = -(m_cluster_sizeX / 2); - ic < (m_cluster_sizeX / 2) + 1; ic++) { - if (ix + ic >= 0 && ix + ic < frame.shape(1) && - iy + ir >= 0 && iy + ir < frame.shape(0)) { - auto tmp = frame(iy + ir, ix + ic) - - pedestal.mean(iy + ir, ix + ic); - cluster.set(i, tmp); - i++; - } - } - } - clusters.push_back(cluster); - } - } - } - return clusters; - } + // if (m_threshold > 2 * tthr) + // tthr = m_threshold - tthr; + // if (m_threshold > 2 * tthr1) + // tthr1 = tthr - tthr1; + // if (m_threshold > 2 * tthr2) + // tthr2 = tthr - tthr2; + // } + // if (total > tthr1 || max > tthr) { + // eventMask[iy][ix] = eventType::PHOTON; + // nph(iy, ix) += 1; + // rest(iy, ix) -= m_threshold; + // } else { + // pedestal.push(iy, ix, frame(iy, ix)); + // continue; + // } + // if (eventMask[iy][ix] == eventType::PHOTON && + // frame(iy, ix) - pedestal.mean(iy, ix) >= max) { + // eventMask[iy][ix] = eventType::PHOTON_MAX; + // DynamicCluster cluster(m_cluster_sizeX, m_cluster_sizeY, + // Dtype(typeid(FRAME_TYPE))); + // cluster.x = ix; + // cluster.y = iy; + // short i = 0; + // for (short ir = -(m_cluster_sizeY / 2); + // ir < (m_cluster_sizeY / 2) + 1; ir++) { + // for (short ic = -(m_cluster_sizeX / 2); + // ic < (m_cluster_sizeX / 2) + 1; ic++) { + // if (ix + ic >= 0 && ix + ic < frame.shape(1) && + // iy + ir >= 0 && iy + ir < frame.shape(0)) { + // auto tmp = frame(iy + ir, ix + ic) - + // pedestal.mean(iy + ir, ix + ic); + // cluster.set(i, tmp); + // i++; + // } + // } + // } + // clusters.push_back(cluster); + // } + // } + // } + // return clusters; + // } }; } // namespace aare \ No newline at end of file diff --git a/include/aare/ClusterVector.hpp b/include/aare/ClusterVector.hpp index 76e7e21..53aeed7 100644 --- a/include/aare/ClusterVector.hpp +++ b/include/aare/ClusterVector.hpp @@ -2,9 +2,18 @@ #include #include #include +#include #include +namespace aare { + +/** + * @brief ClusterVector is a container for clusters of various sizes. It uses a + * contiguous memory buffer to store the clusters. + * @note push_back can invalidate pointers to elements in the container + * @tparam T data type of the pixels in the cluster + */ template class ClusterVector { using value_type = T; using coord_t = int16_t; @@ -24,6 +33,12 @@ template class ClusterVector { constexpr static char m_fmt_base[] = "=h:x:\nh:y:\n({},{}){}:data:" ; public: + /** + * @brief Construct a new ClusterVector object + * @param cluster_size_x size of the cluster in x direction + * @param cluster_size_y size of the cluster in y direction + * @param capacity initial capacity of the buffer in number of clusters + */ ClusterVector(coord_t cluster_size_x, coord_t cluster_size_y, size_t capacity = 1024) : m_cluster_size_x(cluster_size_x), m_cluster_size_y(cluster_size_y), @@ -31,8 +46,6 @@ template class ClusterVector { allocate_buffer(capacity); } ~ClusterVector() { - fmt::print("~ClusterVector - size: {}, capacity: {}\n", m_size, - m_capacity); delete[] m_data; } @@ -63,13 +76,24 @@ template class ClusterVector { return *this; } + /** + * @brief Reserve space for at least capacity clusters + * @param capacity number of clusters to reserve space for + * @note If capacity is less than the current capacity, the function does nothing. + */ void reserve(size_t capacity) { if (capacity > m_capacity) { allocate_buffer(capacity); } } - // data better hold data of the right size! + /** + * @brief Add a cluster to the vector + * @param x x-coordinate of the cluster + * @param y y-coordinate of the cluster + * @param data pointer to the data of the cluster + * @warning The data pointer must point to a buffer of size cluster_size_x * cluster_size_y * sizeof(T) + */ void push_back(coord_t x, coord_t y, const std::byte *data) { if (m_size == m_capacity) { allocate_buffer(m_capacity * 2); @@ -85,6 +109,10 @@ template class ClusterVector { m_size++; } + /** + * @brief Sum the pixels in each cluster + * @return std::vector vector of sums for each cluster + */ std::vector sum() { std::vector sums(m_size); const size_t stride = element_offset(); @@ -101,12 +129,23 @@ template class ClusterVector { } size_t size() const { return m_size; } + size_t capacity() const { return m_capacity; } + + /** + * @brief Return the offset in bytes for a single cluster + */ size_t element_offset() const { return sizeof(m_cluster_size_x) + sizeof(m_cluster_size_y) + m_cluster_size_x * m_cluster_size_y * sizeof(T); } + /** + * @brief Return the offset in bytes for the i-th cluster + */ size_t element_offset(size_t i) const { return element_offset() * i; } + /** + * @brief Return a pointer to the i-th cluster + */ std::byte *element_ptr(size_t i) { return m_data + element_offset(i); } int16_t cluster_size_x() const { return m_cluster_size_x; } @@ -123,13 +162,12 @@ template class ClusterVector { 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 +}; + +} // namespace aare \ No newline at end of file diff --git a/python/src/cluster.hpp b/python/src/cluster.hpp index 0fef093..0467c98 100644 --- a/python/src/cluster.hpp +++ b/python/src/cluster.hpp @@ -10,7 +10,7 @@ #include namespace py = pybind11; -using pd_type = double; +using pd_type = float; template void define_cluster_vector(py::module &m, const std::string &typestr) { @@ -46,7 +46,9 @@ void define_cluster_vector(py::module &m, const std::string &typestr) { void define_cluster_finder_bindings(py::module &m) { py::class_>(m, "ClusterFinder") - .def(py::init, Shape<2>>()) + .def(py::init, Shape<2>, pd_type, size_t>(), py::arg("image_size"), + py::arg("cluster_size"), py::arg("n_sigma") = 5.0, + py::arg("capacity") = 1'000'000) .def("push_pedestal_frame", [](ClusterFinder &self, py::array_t frame) { @@ -66,10 +68,10 @@ void define_cluster_finder_bindings(py::module &m) { return return_image_data(arr); }) .def("steal_clusters", - [](ClusterFinder &self) { - auto v = new ClusterVector(self.steal_clusters()); + [](ClusterFinder &self, bool realloc_same_capacity) { + auto v = new ClusterVector(self.steal_clusters(realloc_same_capacity)); return v; - }) + }, py::arg("realloc_same_capacity") = false) .def("find_clusters", [](ClusterFinder &self, py::array_t frame) { @@ -78,36 +80,11 @@ void define_cluster_finder_bindings(py::module &m) { return; }); - m.def("hello", []() { - fmt::print("Hello from C++\n"); - auto v = new ClusterVector(3, 3); - int data[] = {1, 2, 3, 4, 5, 6, 7, 8, 9}; - v->push_back(5, 30, reinterpret_cast(data)); - v->push_back(5, 55, reinterpret_cast(data)); - v->push_back(5, 20, reinterpret_cast(data)); - v->push_back(5, 30, reinterpret_cast(data)); - - return v; - }); - + define_cluster_vector(m, "i"); define_cluster_vector(m, "d"); + define_cluster_vector(m, "f"); - // py::class_>(m, "ClusterVector", py::buffer_protocol()) - // .def(py::init()) - // .def("size", &ClusterVector::size) - // .def("element_offset", - // py::overload_cast<>(&ClusterVector::element_offset, py::const_)) - // .def_buffer([](ClusterVector &v) -> py::buffer_info { - // return py::buffer_info( - // v.data(), /* Pointer to buffer */ - // v.element_offset(), /* Size of one scalar */ - // fmt::format("h:x:\nh:y:\n({},{})i:data:", v.cluster_size_x(), - // v.cluster_size_y()), /* Format descriptor */ 1, /* Number of - // dimensions */ {v.size()}, /* Buffer dimensions */ - // {v.element_offset()} /* Strides (in bytes) for each index */ - // ); - // }); py::class_(m, "DynamicCluster", py::buffer_protocol()) .def(py::init()) diff --git a/src/ClusterVector.test.cpp b/src/ClusterVector.test.cpp new file mode 100644 index 0000000..ef4e0ee --- /dev/null +++ b/src/ClusterVector.test.cpp @@ -0,0 +1,77 @@ +#include +#include "aare/ClusterVector.hpp" + +// #include +#include + +using aare::ClusterVector; + +TEST_CASE("ClusterVector 2x2 int32_t capacity 4, push back then read") { + struct Cluster_i2x2 { + int16_t x; + int16_t y; + int32_t data[4]; + }; + + ClusterVector cv(2, 2, 4); + REQUIRE(cv.capacity() == 4); + REQUIRE(cv.size() == 0); + REQUIRE(cv.cluster_size_x() == 2); + REQUIRE(cv.cluster_size_y() == 2); + // int16_t, int16_t, 2x2 int32_t = 20 bytes + REQUIRE(cv.element_offset() == 20); + + //Create a cluster and push back into the vector + Cluster_i2x2 c1 = {1, 2, {3, 4, 5, 6}}; + cv.push_back(c1.x, c1.y, reinterpret_cast(&c1.data[0])); + REQUIRE(cv.size() == 1); + REQUIRE(cv.capacity() == 4); + + //Read the cluster back out using copy. TODO! Can we improve the API? + Cluster_i2x2 c2; + std::byte *ptr = cv.element_ptr(0); + std::copy(ptr, ptr + cv.element_offset(), reinterpret_cast(&c2)); + + //Check that the data is the same + REQUIRE(c1.x == c2.x); + REQUIRE(c1.y == c2.y); + for(size_t i = 0; i < 4; i++) { + REQUIRE(c1.data[i] == c2.data[i]); + } +} + +TEST_CASE("Summing 3x1 clusters of int64"){ + struct Cluster_l3x1{ + int16_t x; + int16_t y; + int64_t data[3]; + }; + + ClusterVector cv(3, 1, 2); + REQUIRE(cv.capacity() == 2); + REQUIRE(cv.size() == 0); + REQUIRE(cv.cluster_size_x() == 3); + REQUIRE(cv.cluster_size_y() == 1); + + //Create a cluster and push back into the vector + Cluster_l3x1 c1 = {1, 2, {3, 4, 5}}; + cv.push_back(c1.x, c1.y, reinterpret_cast(&c1.data[0])); + REQUIRE(cv.capacity() == 2); + REQUIRE(cv.size() == 1); + + Cluster_l3x1 c2 = {6, 7, {8, 9, 10}}; + cv.push_back(c2.x, c2.y, reinterpret_cast(&c2.data[0])); + REQUIRE(cv.capacity() == 2); + REQUIRE(cv.size() == 2); + + Cluster_l3x1 c3 = {11, 12, {13, 14, 15}}; + cv.push_back(c3.x, c3.y, reinterpret_cast(&c3.data[0])); + REQUIRE(cv.capacity() == 4); + REQUIRE(cv.size() == 3); + + auto sums = cv.sum(); + REQUIRE(sums.size() == 3); + REQUIRE(sums[0] == 12); + REQUIRE(sums[1] == 27); + REQUIRE(sums[2] == 42); +} \ No newline at end of file