diff --git a/CMakeLists.txt b/CMakeLists.txt index 9966bb4..dbda663 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -55,6 +55,10 @@ option( "Install the python extension in the install tree under CMAKE_INSTALL_PREFIX/aare/" OFF) option(AARE_ASAN "Enable AddressSanitizer" OFF) +option( + AARE_TUNE_LOCAL + "Optimize for the building machine's CPU (-march=native -mtune=native). Not portable, the resulting binaries might not run on other machines." + OFF) # Configure which of the dependencies to use FetchContent for option(AARE_FETCH_FMT "Use FetchContent to download fmt" ON) @@ -352,6 +356,22 @@ else() target_compile_options(aare_compiler_flags INTERFACE -Werror) endif() + if(AARE_TUNE_LOCAL) + include(CheckCXXCompilerFlag) + check_cxx_compiler_flag("-march=native" AARE_HAS_MARCH_NATIVE) + check_cxx_compiler_flag("-mtune=native" AARE_HAS_MTUNE_NATIVE) + if(AARE_HAS_MARCH_NATIVE AND AARE_HAS_MTUNE_NATIVE) + message(STATUS "Tuning for the local CPU: -march=native -mtune=native") + target_compile_options(aare_compiler_flags INTERFACE -march=native + -mtune=native) + else() + message( + WARNING + "AARE_TUNE_LOCAL requested but the compiler does not support -march=native/-mtune=native. Ignoring." + ) + endif() + endif() + endif() # GCC/Clang specific if(AARE_PYTHON_BINDINGS) diff --git a/include/aare/ClusterFinder.hpp b/include/aare/ClusterFinder.hpp index c3089b7..af41cda 100644 --- a/include/aare/ClusterFinder.hpp +++ b/include/aare/ClusterFinder.hpp @@ -74,11 +74,7 @@ class ClusterFinder { NDArray noise() { return m_pedestal.std(); } void clear_pedestal() { m_pedestal.clear(); } - /** - * @brief Refresh the cached std of the underlying pedestal. Call before - * reading the pedestal's cached std. - */ - void update_std() { m_pedestal.update_std(); } + void update_threshold() { m_threshold = m_pedestal.std() * m_nSigma; } @@ -161,8 +157,7 @@ class ClusterFinder { } else if (total > c3 * threshold) { // pass, store the cluster below } else { - // m_pedestal.push(iy, ix, frame(iy, ix)); // Safe option - m_pedestal.push(iy, ix, frame.data()[center]); + m_pedestal.push_fast(center, frame.data()[center]); return; // It was a pedestal value nothing to store } @@ -254,7 +249,7 @@ class ClusterFinder { auto pd = m_pedestal.view().data(); auto corrected = m_pd_corrected_frame.data(); auto frame_data = frame.data(); - for (size_t i = 0; i < n_pixels; i++) { + for (ssize_t i = 0; i < n_pixels; i++) { corrected[i] = static_cast(frame_data[i]) - pd[i]; } diff --git a/include/aare/FastPedestal.hpp b/include/aare/FastPedestal.hpp index 1c17c48..b371be9 100644 --- a/include/aare/FastPedestal.hpp +++ b/include/aare/FastPedestal.hpp @@ -4,6 +4,7 @@ #include "aare/NDArray.hpp" #include "aare/NDView.hpp" #include +#include namespace aare { @@ -14,10 +15,9 @@ namespace aare { * * @tparam PEDESTAL_TYPE type of the exposed mean and std */ -template class FastPedestal { - // TODO! Force floating point sum type? - // how does the internal calculation work with integers? +template class FastPedestal { + // Did we accumulate enough samples and updated the mean? bool m_ready = false; uint32_t m_rows; @@ -25,62 +25,71 @@ template class FastPedestal { uint32_t m_samples; double m_inv_samples; // precompute 1/m_samples for faster division - uint32_t m_cur_samples = 0; // TODO! do we need this when we have m_samples? + uint32_t m_cur_samples = 0; // number of samples accumulated so far - // for cache we want to keep sum and sum2 close + // For cache locality we want to keep sum and sum2 close. Improves performance + // for random access. struct Entry { double sum; double sum2; }; - // TODO! in case of int needs to be changed to uint64_t NDArray m_sum; // Cache mean since it is used over and over in the ClusterFinder // This optimization is related to the access pattern of the ClusterFinder // Relies on having more reads than pushes to the pedestal + // But also makes sense when subtracting the pedestal from the frame NDArray m_mean; - // Cache std. Only refreshed via update_std() to keep push() cheap. - // This gives a measureable speedup in the ClusterFinder. - NDArray m_std; - + // Helper function to convert row and column indices to a flat index + // used to provide both row column and flat index access to the pedestal + size_t rc_to_index(uint32_t row, uint32_t col) const { + return (static_cast(row) * m_cols) + col; + } public: FastPedestal(uint32_t rows, uint32_t cols, uint32_t n_samples = 1000) : m_rows(rows), m_cols(cols), m_samples(n_samples), m_inv_samples(1.0 / n_samples), m_sum({rows, cols}, Entry{0, 0}), - m_mean({rows, cols}), - m_std({rows, cols}, 0) { + m_mean({rows, cols}, PEDESTAL_TYPE(0)) { assert(rows > 0 && cols > 0 && n_samples > 0); } ~FastPedestal() = default; - NDArray mean() { return m_mean; } - const NDView view() const { return m_mean.view(); } - PEDESTAL_TYPE mean(const uint32_t row, const uint32_t col) const { + NDArray mean() { return m_mean; } + + PEDESTAL_TYPE mean(uint32_t row, uint32_t col) const { return m_mean(row, col); } - NDArray cached_std() { return m_std; } - - PEDESTAL_TYPE cached_std(const uint32_t row, const uint32_t col) const { - return m_std(row, col); - } - - PEDESTAL_TYPE variance(const uint32_t row, const uint32_t col) const { - auto &entry = m_sum(row, col); - auto m2 = entry.sum * m_inv_samples * entry.sum * m_inv_samples; - return entry.sum2 * m_inv_samples - m2; + PEDESTAL_TYPE mean(ssize_t index) const { + return m_mean[index]; } NDArray variance() { NDArray res({m_rows, m_cols}); - for (ssize_t row = 0; row < m_rows; ++row) { - for (ssize_t col = 0; col < m_cols; ++col) { - res(row, col) = variance(row, col); - } + for (ssize_t i = 0; i < m_sum.size(); ++i) { + res[i] = variance(i); + } + return res; + } + + PEDESTAL_TYPE variance(const uint32_t row, const uint32_t col) const { + return variance(rc_to_index(row, col)); + } + + PEDESTAL_TYPE variance(ssize_t index) const { + auto &entry = m_sum[index]; + auto m2 = entry.sum * m_inv_samples * entry.sum * m_inv_samples; + return entry.sum2 * m_inv_samples - m2; + } + + NDArray std() { + NDArray res({m_rows, m_cols}); + for (ssize_t i = 0; i < m_sum.size(); ++i) { + res[i] = std(i); } return res; } @@ -89,32 +98,23 @@ template class FastPedestal { return std::sqrt(variance(row, col)); } - NDArray std() { - NDArray res({m_rows, m_cols}); - for (ssize_t row = 0; row < m_rows; ++row) { - for (ssize_t col = 0; col < m_cols; ++col) { - res(row, col) = std(row, col); - } - } - return res; + PEDESTAL_TYPE std(ssize_t index) const { + return std::sqrt(variance(index)); } - bool ready() { return m_ready; } - uint32_t cur_samples() { return m_cur_samples; } + bool ready() const { return m_ready; } + + uint32_t cur_samples() const { return m_cur_samples; } void clear() { - m_sum = Entry{double(0), double(0)}; - m_mean = PEDESTAL_TYPE(0); - m_std = PEDESTAL_TYPE(0); + m_sum = Entry{0., 0.}; + m_mean = PEDESTAL_TYPE(0.); m_ready = false; } template void push(NDView frame) { - assert(frame.size() == m_rows * m_cols); - - // TODO! move away from m_rows, m_cols if (frame.shape() != std::array{m_rows, m_cols}) { throw std::runtime_error( "Frame shape does not match pedestal shape"); @@ -127,9 +127,6 @@ template class FastPedestal { } } template void push_init(NDView frame) { - assert(frame.size() == m_rows * m_cols); - - // TODO! move away from m_rows, m_cols if (frame.shape() != std::array{m_rows, m_cols}) { throw std::runtime_error( "Frame shape does not match pedestal shape"); @@ -151,7 +148,6 @@ template class FastPedestal { m_cur_samples += 1; if (m_cur_samples == m_samples) { - update_std(); update_mean(); m_ready = true; } @@ -168,18 +164,35 @@ template class FastPedestal { uint32_t cols() const { return m_cols; } uint32_t n_samples() const { return m_samples; } - // pixel level operations (should be refactored to allow users to implement - // their own pixel level operations) + + /** + * @brief Update one pixel using its flat index. + * + * This steady-state fast path assumes the pedestal is ready and the index + * is valid. Assertions check those preconditions in debug builds. + */ + template + void push_fast(const std::size_t index, const T value) noexcept { + assert(m_ready); + assert(index < static_cast(m_sum.size())); + + const auto val = static_cast(value); + auto &entry = m_sum[index]; + entry.sum += val - entry.sum * m_inv_samples; + entry.sum2 += val * val - entry.sum2 * m_inv_samples; + m_mean[index] = + static_cast(entry.sum * m_inv_samples); + } + + template void push(const uint32_t row, const uint32_t col, const T val_) { if (!ready()) { throw std::runtime_error("Pedestal is not ready, cannot push"); } - auto val = static_cast(val_); - auto &entry = m_sum(row, col); - entry.sum += val - entry.sum * m_inv_samples; - entry.sum2 += val * val - entry.sum2 * m_inv_samples; - m_mean(row, col) = entry.sum * m_inv_samples; + const auto index = + (static_cast(row) * m_cols) + col; + push_fast(index, val_); } template @@ -198,25 +211,11 @@ template class FastPedestal { * push_no_update. It is not necessary to call this function after push. */ void update_mean() { - for (size_t row = 0; row < m_rows; row++) { - for (size_t col = 0; col < m_cols; col++) { - const auto &entry = m_sum(row, col); - m_mean(row, col) = entry.sum * m_inv_samples; - } + for (ssize_t i = 0; i < m_sum.size(); i++) { + auto &entry = m_sum[i]; + m_mean[i] = static_cast(entry.sum * m_inv_samples); } } - /** - * @brief Refresh the cached std for all pixels from the current sums. - * Kept separate from push() so that pushes stay cheap; call this before - * reading cached_std(). - */ - void update_std() { - for (size_t row = 0; row < m_rows; row++) { - for (size_t col = 0; col < m_cols; col++) { - m_std(row, col) = std(row, col); - } - } - } }; } // namespace aare diff --git a/python/src/bind_ClusterFinder.hpp b/python/src/bind_ClusterFinder.hpp index 564e0cf..fa97915 100644 --- a/python/src/bind_ClusterFinder.hpp +++ b/python/src/bind_ClusterFinder.hpp @@ -50,8 +50,6 @@ void define_ClusterFinder(py::module &m, const std::string &typestr) { }) .def("clear_pedestal", &ClusterFinder::clear_pedestal) - .def("update_std", - &ClusterFinder::update_std) .def("update_threshold", &ClusterFinder::update_threshold) .def_property_readonly( diff --git a/python/src/fast_pedestal.hpp b/python/src/fast_pedestal.hpp index 73651b3..7596455 100644 --- a/python/src/fast_pedestal.hpp +++ b/python/src/fast_pedestal.hpp @@ -51,12 +51,6 @@ void define_fast_pedestal_bindings(py::module &m, const std::string &name) { *standard_deviation = self.std(); return return_image_data(standard_deviation); }) - .def("cached_std", - [](FastPedestal &self) { - auto standard_deviation = new NDArray{}; - *standard_deviation = self.cached_std(); - return return_image_data(standard_deviation); - }) .def( "__array_ufunc__", [](py::object self, py::object ufunc, const std::string &method, @@ -100,15 +94,8 @@ void define_fast_pedestal_bindings(py::module &m, const std::string &name) { pedestal.push_init(make_view_2d(frame)); }, py::arg("frame").noconvert()) - // .def( - // "push_no_update", - // [](FastPedestal &pedestal, - // py::array_t &frame) { - // pedestal.push_no_update(make_view_2d(frame)); - // }, - // py::arg("frame").noconvert()) + .def("update_mean", &FastPedestal::update_mean) - .def("update_std", &FastPedestal::update_std) .def_buffer([](FastPedestal &self) { auto mean = self.view(); return py::buffer_info( diff --git a/python/src/pedestal.hpp b/python/src/pedestal.hpp index e38ba62..229f3ad 100644 --- a/python/src/pedestal.hpp +++ b/python/src/pedestal.hpp @@ -110,7 +110,6 @@ void define_pedestal_bindings(py::module &m, const std::string &name) { }, py::arg().noconvert()) .def("update_mean", &Pedestal::update_mean) - .def("update_std", &Pedestal::update_std) .def_buffer([](Pedestal &self) { auto mean = self.view(); return py::buffer_info(