From 13471582354bdf235bd65188f00c6f9cf19e0f67 Mon Sep 17 00:00:00 2001 From: Alice Date: Thu, 24 Jul 2025 15:40:05 +0200 Subject: [PATCH 1/5] templated calculate_pedestal with boolean template argument only_gain0, added drop_dimension to NDArray and reference pointer to data --- CMakeLists.txt | 1 + include/aare/NDArray.hpp | 33 +++++++++---- include/aare/calibration.hpp | 83 +++++++++++++++++---------------- python/src/bind_calibration.hpp | 4 +- src/calibration.cpp | 56 ++++++---------------- 5 files changed, 85 insertions(+), 92 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 3cfc6a7..6338f93 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -438,6 +438,7 @@ endif() if(AARE_TESTS) set(TestSources ${CMAKE_CURRENT_SOURCE_DIR}/src/algorithm.test.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/src/calibration.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/defs.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/decode.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/Dtype.test.cpp diff --git a/include/aare/NDArray.hpp b/include/aare/NDArray.hpp index f73bfd6..c3edbfb 100644 --- a/include/aare/NDArray.hpp +++ b/include/aare/NDArray.hpp @@ -33,7 +33,7 @@ class NDArray : public ArrayExpr, Ndim> { * @brief Default constructor. Will construct an empty NDArray. * */ - NDArray() : shape_(), strides_(c_strides(shape_)), data_(nullptr){}; + NDArray() : shape_(), strides_(c_strides(shape_)), data_(nullptr) {}; /** * @brief Construct a new NDArray object with a given shape. @@ -185,6 +185,7 @@ class NDArray : public ArrayExpr, Ndim> { const T &operator[](ssize_t i) const { return data_[i]; } T *data() { return data_; } + T *&data_ref() { return data_; } std::byte *buffer() { return reinterpret_cast(data_); } ssize_t size() const { return static_cast(size_); } size_t total_bytes() const { return size_ * sizeof(T); } @@ -211,13 +212,30 @@ class NDArray : public ArrayExpr, Ndim> { void Print_all(); void Print_some(); + template 2)>> + NDArray drop_dimension() && { + + std::array new_shape; + + std::copy(shape_.begin() + 1, shape_.begin() + Ndim, new_shape.begin()); + + NDArray new_array(new_shape); + + delete new_array.data(); + + new_array.data_ref() = data_; + + this->reset(); + + return new_array; + } + void reset() { data_ = nullptr; size_ = 0; std::fill(shape_.begin(), shape_.end(), 0); std::fill(strides_.begin(), strides_.end(), 0); } - }; // Move assign @@ -382,8 +400,6 @@ NDArray NDArray::operator*(const T &value) { return result; } - - template std::ostream &operator<<(std::ostream &os, const NDArray &arr) { for (auto row = 0; row < arr.shape(0); ++row) { @@ -434,17 +450,18 @@ NDArray load(const std::string &pathname, return img; } - template NDArray safe_divide(const NDArray &numerator, - const NDArray &denominator) { + const NDArray &denominator) { if (numerator.shape() != denominator.shape()) { - throw std::runtime_error("Shapes of numerator and denominator must match"); + throw std::runtime_error( + "Shapes of numerator and denominator must match"); } NDArray result(numerator.shape()); for (ssize_t i = 0; i < numerator.size(); ++i) { if (denominator[i] != 0) { - result[i] = static_cast(numerator[i]) / static_cast(denominator[i]); + result[i] = + static_cast(numerator[i]) / static_cast(denominator[i]); } else { result[i] = RT{0}; // or handle division by zero as needed } diff --git a/include/aare/calibration.hpp b/include/aare/calibration.hpp index bf543aa..a94705f 100644 --- a/include/aare/calibration.hpp +++ b/include/aare/calibration.hpp @@ -103,7 +103,7 @@ void apply_calibration_impl(NDView res, NDView raw_data, } } -template +template void apply_calibration(NDView res, NDView raw_data, NDView ped, NDView cal, ssize_t n_threads = 4) { @@ -113,52 +113,44 @@ void apply_calibration(NDView res, NDView raw_data, for (const auto &lim : limits) futures.push_back(std::async( static_cast, NDView, - NDView, NDView, int, - int)>( + NDView, NDView, int, int)>( apply_calibration_impl), res, raw_data, ped, cal, lim.first, lim.second)); for (auto &f : futures) f.get(); } - +template std::pair, NDArray> -sum_and_count_per_gain(NDView raw_data); - -std::pair, NDArray> -sum_and_count_g0(NDView raw_data); - -template -NDArray calculate_pedestal_g0(NDView raw_data, - ssize_t n_threads) { - std::vector, NDArray>>> - futures; - futures.reserve(n_threads); - - auto subviews = make_subviews(raw_data, n_threads); - - for (auto view : subviews) { - futures.push_back(std::async( - static_cast, NDArray> (*)( - NDView)>(&sum_and_count_g0), - view)); - } - NDArray accumulator( - std::array{raw_data.shape(1), raw_data.shape(2)}, 0); - NDArray count( - std::array{raw_data.shape(1), raw_data.shape(2)}, 0); - for (auto &f : futures) { - auto [acc, cnt] = f.get(); - accumulator += acc; - count += cnt; +sum_and_count_per_gain(NDView raw_data) { + constexpr ssize_t num_gains = only_gain0 ? 1 : 3; + NDArray accumulator( + std::array{num_gains, raw_data.shape(1), raw_data.shape(2)}, + 0); + NDArray count( + std::array{num_gains, raw_data.shape(1), raw_data.shape(2)}, + 0); + for (int frame_nr = 0; frame_nr != raw_data.shape(0); ++frame_nr) { + for (int row = 0; row != raw_data.shape(1); ++row) { + for (int col = 0; col != raw_data.shape(2); ++col) { + auto [value, gain] = + get_value_and_gain(raw_data(frame_nr, row, col)); + if (gain != 0 && only_gain0) + continue; + accumulator(gain, row, col) += value; + count(gain, row, col) += 1; + } + } } - return safe_divide(accumulator, count); + return {std::move(accumulator), std::move(count)}; } -template -NDArray calculate_pedestal(NDView raw_data, - ssize_t n_threads) { +template +NDArray(only_gain0)> +calculate_pedestal(NDView raw_data, ssize_t n_threads) { + + constexpr ssize_t num_gains = only_gain0 ? 1 : 3; std::vector, NDArray>>> futures; futures.reserve(n_threads); @@ -168,21 +160,27 @@ NDArray calculate_pedestal(NDView raw_data, for (auto view : subviews) { futures.push_back(std::async( static_cast, NDArray> (*)( - NDView)>(&sum_and_count_per_gain), + NDView)>(&sum_and_count_per_gain), view)); } NDArray accumulator( - std::array{3, raw_data.shape(1), raw_data.shape(2)}, 0); + std::array{num_gains, raw_data.shape(1), raw_data.shape(2)}, + 0); NDArray count( - std::array{3, raw_data.shape(1), raw_data.shape(2)}, 0); + std::array{num_gains, raw_data.shape(1), raw_data.shape(2)}, + 0); for (auto &f : futures) { auto [acc, cnt] = f.get(); accumulator += acc; count += cnt; } - return safe_divide(accumulator, count); + if constexpr (only_gain0) { + return safe_divide(accumulator, count).drop_dimension(); + } else { + return safe_divide(accumulator, count); + } } /** @@ -205,4 +203,9 @@ NDArray count_switching_pixels(NDView raw_data); NDArray count_switching_pixels(NDView raw_data, ssize_t n_threads); +template +auto calculate_pedestal_g0(NDView raw_data, ssize_t n_threads) { + return calculate_pedestal(raw_data, n_threads); +} + } // namespace aare \ No newline at end of file diff --git a/python/src/bind_calibration.hpp b/python/src/bind_calibration.hpp index b95fa34..836a6de 100644 --- a/python/src/bind_calibration.hpp +++ b/python/src/bind_calibration.hpp @@ -56,7 +56,7 @@ py::array_t pybind_calculate_pedestal( auto data_span = make_view_3d(data); auto arr = new NDArray{}; - *arr = aare::calculate_pedestal(data_span, n_threads); + *arr = aare::calculate_pedestal(data_span, n_threads); return return_image_data(arr); } @@ -67,7 +67,7 @@ py::array_t pybind_calculate_pedestal_g0( auto data_span = make_view_3d(data); auto arr = new NDArray{}; - *arr = aare::calculate_pedestal_g0(data_span, n_threads); + *arr = aare::calculate_pedestal(data_span, n_threads); return return_image_data(arr); } diff --git a/src/calibration.cpp b/src/calibration.cpp index eeaae84..0cb3b99 100644 --- a/src/calibration.cpp +++ b/src/calibration.cpp @@ -2,46 +2,14 @@ namespace aare { -std::pair, NDArray> sum_and_count_per_gain(NDView raw_data){ - NDArray accumulator(std::array{3, raw_data.shape(1), raw_data.shape(2)}, 0); - NDArray count(std::array{3, raw_data.shape(1), raw_data.shape(2)}, 0); - for (int frame_nr = 0; frame_nr != raw_data.shape(0); ++frame_nr) { - for (int row = 0; row != raw_data.shape(1); ++row) { - for (int col = 0; col != raw_data.shape(2); ++col) { - auto [value, gain] = get_value_and_gain(raw_data(frame_nr, row, col)); - accumulator(gain, row, col) += value; - count(gain, row, col) += 1; - } - } - } - - return {std::move(accumulator), std::move(count)}; -} - -std::pair, NDArray> sum_and_count_g0(NDView raw_data){ - NDArray accumulator(std::array{raw_data.shape(1), raw_data.shape(2)}, 0); - NDArray count(std::array{raw_data.shape(1), raw_data.shape(2)}, 0); - for (int frame_nr = 0; frame_nr != raw_data.shape(0); ++frame_nr) { - for (int row = 0; row != raw_data.shape(1); ++row) { - for (int col = 0; col != raw_data.shape(2); ++col) { - auto [value, gain] = get_value_and_gain(raw_data(frame_nr, row, col)); - if (gain != 0) - continue; // we only care about gain 0 - accumulator(row, col) += value; - count(row, col) += 1; - } - } - } - - return {std::move(accumulator), std::move(count)}; -} - NDArray count_switching_pixels(NDView raw_data) { - NDArray switched(std::array{raw_data.shape(1), raw_data.shape(2)},0); + NDArray switched( + std::array{raw_data.shape(1), raw_data.shape(2)}, 0); for (int frame_nr = 0; frame_nr != raw_data.shape(0); ++frame_nr) { for (int row = 0; row != raw_data.shape(1); ++row) { for (int col = 0; col != raw_data.shape(2); ++col) { - auto [value, gain] = get_value_and_gain(raw_data(frame_nr, row, col)); + auto [value, gain] = + get_value_and_gain(raw_data(frame_nr, row, col)); if (gain != 0) { switched(row, col) += 1; } @@ -51,22 +19,26 @@ NDArray count_switching_pixels(NDView raw_data) { return switched; } - -NDArray count_switching_pixels(NDView raw_data, ssize_t n_threads){ - NDArray switched(std::array{raw_data.shape(1), raw_data.shape(2)},0); +NDArray count_switching_pixels(NDView raw_data, + ssize_t n_threads) { + NDArray switched( + std::array{raw_data.shape(1), raw_data.shape(2)}, 0); std::vector>> futures; futures.reserve(n_threads); auto subviews = make_subviews(raw_data, n_threads); for (auto view : subviews) { - futures.push_back(std::async(static_cast(*)(NDView)>(&count_switching_pixels), view)); + futures.push_back( + std::async(static_cast (*)(NDView)>( + &count_switching_pixels), + view)); } - + for (auto &f : futures) { switched += f.get(); } return switched; } -}// namespace aare \ No newline at end of file +} // namespace aare \ No newline at end of file From 1195a5e100b355c9b6353e2bc0c5299547e6ab58 Mon Sep 17 00:00:00 2001 From: Alice Date: Fri, 25 Jul 2025 10:18:55 +0200 Subject: [PATCH 2/5] added drop dimension test, added file calibration.test.cpp --- src/NDArray.test.cpp | 19 +++++++++++++++- src/calibration.test.cpp | 49 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 67 insertions(+), 1 deletion(-) create mode 100644 src/calibration.test.cpp diff --git a/src/NDArray.test.cpp b/src/NDArray.test.cpp index 9af8379..97a7f09 100644 --- a/src/NDArray.test.cpp +++ b/src/NDArray.test.cpp @@ -427,4 +427,21 @@ TEST_CASE("Construct an NDArray from an std::array") { for (uint32_t i = 0; i < a.size(); ++i) { REQUIRE(a(i) == b[i]); } -} \ No newline at end of file +} + +TEST_CASE("Drop dimension") { + NDArray array(std::array{2, 2, 2}); + + std::fill(array.begin(), array.begin() + 4, 0); + std::fill(array.begin() + 4, array.end(), 1); + + auto new_array = std::move(array).drop_dimension(); + + CHECK(new_array.shape() == std::array{2, 2}); + CHECK(new_array.size() == 4); + + std::for_each(new_array.begin(), new_array.end(), + [](int &element) { CHECK(element == 0); }); + + CHECK(array.size() == 0); // array was moved +} diff --git a/src/calibration.test.cpp b/src/calibration.test.cpp new file mode 100644 index 0000000..712dcd8 --- /dev/null +++ b/src/calibration.test.cpp @@ -0,0 +1,49 @@ +/************************************************ + * @file test-Cluster.cpp + * @short test case for generic Cluster, ClusterVector, and calculate_eta2 + ***********************************************/ + +#include "aare/calibration.hpp" + +// #include "catch.hpp" +#include +#include +#include + +using namespace aare; + +TEST_CASE("Test Pedestal Generation", "[.calibration]") { + NDArray raw(std::array{3, 2, 2}, 0); + + // gain 0 + raw(0, 0, 0) = 100; + raw(1, 0, 0) = 200; + raw(2, 0, 0) = 300; + + // gain 1 + raw(0, 0, 1) = (1 << 14) + 100; + raw(1, 0, 1) = (1 << 14) + 200; + raw(2, 0, 1) = (1 << 14) + 300; + + raw(0, 1, 0) = (1 << 14) + 37; + raw(1, 1, 0) = 38; + raw(2, 1, 0) = (3 << 14) + 39; + + // gain 2 + raw(0, 1, 1) = (3 << 14) + 100; + raw(1, 1, 1) = (3 << 14) + 200; + raw(2, 1, 1) = (3 << 14) + 300; + + auto pedestal = calculate_pedestal(raw.view(), 4); + + REQUIRE(pedestal.size() == raw.size()); + CHECK(pedestal(0, 0, 0) == 200); + CHECK(pedestal(1, 0, 0) == 0); + CHECK(pedestal(1, 0, 1) == 200); + + auto pedestal_gain0 = calculate_pedestal_g0(raw.view(), 4); + + REQUIRE(pedestal_gain0.size() == 4); + CHECK(pedestal_gain0(0, 0) == 200); + CHECK(pedestal_gain0(1, 0) == 38); +} \ No newline at end of file From d6222027d07fe7bc6f89661f856301ec970de633 Mon Sep 17 00:00:00 2001 From: froejdh_e Date: Fri, 25 Jul 2025 10:40:32 +0200 Subject: [PATCH 3/5] move constructor for Ndim-1 --- include/aare/NDArray.hpp | 41 ++++++++++++++++++------------------ include/aare/NDView.hpp | 26 +++++++++++++++++++++++ include/aare/calibration.hpp | 20 ++++++++---------- src/NDArray.test.cpp | 24 +++++++++++++++++++++ 4 files changed, 79 insertions(+), 32 deletions(-) diff --git a/include/aare/NDArray.hpp b/include/aare/NDArray.hpp index c3edbfb..3495276 100644 --- a/include/aare/NDArray.hpp +++ b/include/aare/NDArray.hpp @@ -25,7 +25,7 @@ template class NDArray : public ArrayExpr, Ndim> { std::array shape_; std::array strides_; - size_t size_{}; + size_t size_{}; //TODO! do we need to store size when we have shape? T *data_; public: @@ -43,8 +43,7 @@ class NDArray : public ArrayExpr, Ndim> { */ explicit NDArray(std::array shape) : shape_(shape), strides_(c_strides(shape_)), - size_(std::accumulate(shape_.begin(), shape_.end(), 1, - std::multiplies<>())), + size_(num_elements(shape_)), data_(new T[size_]) {} /** @@ -79,6 +78,24 @@ class NDArray : public ArrayExpr, Ndim> { other.reset(); // TODO! is this necessary? } + + //Move constructor from an an array with Ndim + 1 + template > + NDArray(NDArray &&other) + : shape_(drop_first_dim(other.shape())), + strides_(c_strides(shape_)), size_(num_elements(shape_)), + data_(other.data()) { + + // For now only allow move if the size matches, to avoid unreachable data + // if the use case arises we can remove this check + if(size() != other.size()) { + data_ = nullptr; // avoid double free, other will clean up the memory in it's destructor + throw std::runtime_error(LOCATION + + "Size mismatch in move constructor of NDArray"); + } + other.reset(); + } + // Copy constructor NDArray(const NDArray &other) : shape_(other.shape_), strides_(c_strides(shape_)), @@ -212,24 +229,6 @@ class NDArray : public ArrayExpr, Ndim> { void Print_all(); void Print_some(); - template 2)>> - NDArray drop_dimension() && { - - std::array new_shape; - - std::copy(shape_.begin() + 1, shape_.begin() + Ndim, new_shape.begin()); - - NDArray new_array(new_shape); - - delete new_array.data(); - - new_array.data_ref() = data_; - - this->reset(); - - return new_array; - } - void reset() { data_ = nullptr; size_ = 0; diff --git a/include/aare/NDView.hpp b/include/aare/NDView.hpp index 56bca2f..c946c26 100644 --- a/include/aare/NDView.hpp +++ b/include/aare/NDView.hpp @@ -26,6 +26,32 @@ Shape make_shape(const std::vector &shape) { return arr; } + +/** + * @brief Helper function to drop the first dimension of a shape. + * This is useful when you want to create a 2D view from a 3D array. + * @param shape The shape to drop the first dimension from. + * @return A new shape with the first dimension dropped. + */ +template +Shape drop_first_dim(const Shape &shape) { + Shape new_shape; + std::copy(shape.begin() + 1, shape.end(), new_shape.begin()); + return new_shape; +} + +/** + * @brief Helper function when constructing NDArray/NDView. Calculates the number + * of elements in the resulting array from a shape. + * @param shape The shape to calculate the number of elements for. + * @return The number of elements in and NDArray/NDView of that shape. + */ +template +size_t num_elements(const Shape &shape) { + return std::accumulate(shape.begin(), shape.end(), 1, + std::multiplies()); +} + template ssize_t element_offset(const Strides & /*unused*/) { return 0; diff --git a/include/aare/calibration.hpp b/include/aare/calibration.hpp index a94705f..3b18da6 100644 --- a/include/aare/calibration.hpp +++ b/include/aare/calibration.hpp @@ -163,24 +163,22 @@ calculate_pedestal(NDView raw_data, ssize_t n_threads) { NDView)>(&sum_and_count_per_gain), view)); } + Shape<3> shape{num_gains, raw_data.shape(1), raw_data.shape(2)}; + NDArray accumulator(shape, 0); + NDArray count(shape, 0); - NDArray accumulator( - std::array{num_gains, raw_data.shape(1), raw_data.shape(2)}, - 0); - NDArray count( - std::array{num_gains, raw_data.shape(1), raw_data.shape(2)}, - 0); + // Combine the results from the futures for (auto &f : futures) { auto [acc, cnt] = f.get(); accumulator += acc; count += cnt; } - if constexpr (only_gain0) { - return safe_divide(accumulator, count).drop_dimension(); - } else { - return safe_divide(accumulator, count); - } + + // Will move to a NDArray(only_gain0)> + // if only_gain0 is true + return safe_divide(accumulator, count); + } /** diff --git a/src/NDArray.test.cpp b/src/NDArray.test.cpp index 9af8379..a4d687a 100644 --- a/src/NDArray.test.cpp +++ b/src/NDArray.test.cpp @@ -427,4 +427,28 @@ TEST_CASE("Construct an NDArray from an std::array") { for (uint32_t i = 0; i < a.size(); ++i) { REQUIRE(a(i) == b[i]); } +} + + +TEST_CASE("Move construct from an array with Ndim + 1") { + NDArray a({{1,2,2}}, 0); + a(0, 0, 0) = 1; + a(0, 0, 1) = 2; + a(0, 1, 0) = 3; + a(0, 1, 1) = 4; + + + NDArray b(std::move(a)); + REQUIRE(b.shape() == Shape<2>{2,2}); + REQUIRE(b.size() == 4); + REQUIRE(b(0, 0) == 1); + REQUIRE(b(0, 1) == 2); + REQUIRE(b(1, 0) == 3); + REQUIRE(b(1, 1) == 4); + +} + +TEST_CASE("Move construct from an array with Ndim + 1 throws on size mismatch") { + NDArray a({{2,2,2}}, 0); + REQUIRE_THROWS(NDArray(std::move(a))); } \ No newline at end of file From 3d6858ad3382a6ee11abc6357697382eac7f94f9 Mon Sep 17 00:00:00 2001 From: froejdh_e Date: Fri, 25 Jul 2025 10:42:47 +0200 Subject: [PATCH 4/5] removed data_ref --- include/aare/NDArray.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/include/aare/NDArray.hpp b/include/aare/NDArray.hpp index 3495276..731ba2f 100644 --- a/include/aare/NDArray.hpp +++ b/include/aare/NDArray.hpp @@ -202,7 +202,6 @@ class NDArray : public ArrayExpr, Ndim> { const T &operator[](ssize_t i) const { return data_[i]; } T *data() { return data_; } - T *&data_ref() { return data_; } std::byte *buffer() { return reinterpret_cast(data_); } ssize_t size() const { return static_cast(size_); } size_t total_bytes() const { return size_ * sizeof(T); } From 89bb8776ea8ff1699ec340bf279e05a78b503687 Mon Sep 17 00:00:00 2001 From: froejdh_e Date: Fri, 25 Jul 2025 11:44:27 +0200 Subject: [PATCH 5/5] check Ndim on drop_first_dim --- include/aare/NDView.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/include/aare/NDView.hpp b/include/aare/NDView.hpp index c946c26..0aa4f78 100644 --- a/include/aare/NDView.hpp +++ b/include/aare/NDView.hpp @@ -35,6 +35,7 @@ Shape make_shape(const std::vector &shape) { */ template Shape drop_first_dim(const Shape &shape) { + static_assert(Ndim > 1, "Cannot drop first dimension from a 1D shape"); Shape new_shape; std::copy(shape.begin() + 1, shape.end(), new_shape.begin()); return new_shape;