diff --git a/include/aare/ClusterFinderCUDA.hpp b/include/aare/ClusterFinderCUDA.hpp index 42cd3312..2fdb7c65 100644 --- a/include/aare/ClusterFinderCUDA.hpp +++ b/include/aare/ClusterFinderCUDA.hpp @@ -133,7 +133,7 @@ class ClusterFinderCUDA { size_t m_output_bytes_per_frame; COMPUTE_TYPE m_nSigma; - Pedestal m_pedestal; + FastPedestal m_pedestal; ClusterVector m_clusters; bool m_pedestal_dirty = true; @@ -420,6 +420,10 @@ class ClusterFinderCUDA { * fixed-size D2H * @param n_streams_ number of CUDA streams for multi-frame * overlap + * @param min_pedestal_samples number of frames the host pedestal + * accumulates before it is ready, and the reciprocal of the steady-state + * update weight. Matches ClusterFinder so a CPU/GPU comparison warms up + * identically. * @param time_kernels enable per-frame CUDA-event kernel * timing. Off by default: it adds two event records per frame to the * streams and one host-side query per frame, and the resulting number is @@ -427,12 +431,13 @@ class ClusterFinderCUDA { */ ClusterFinderCUDA(Shape<2> shape_, COMPUTE_TYPE nSigma = 5.0, size_t max_clusters_per_frame = 2048, int n_streams_ = 4, + size_t min_pedestal_samples = 1000, bool time_kernels = false) : m_shape(shape_), nrows(shape_[0]), ncols(shape_[1]), m_image_size(nrows * ncols), n_streams(n_streams_), m_max_clusters_per_frame(max_clusters_per_frame), m_nSigma(nSigma), - m_pedestal(shape_[0], shape_[1]), m_clusters(max_clusters_per_frame), - m_time_kernels(time_kernels) { + m_pedestal(shape_[0], shape_[1], min_pedestal_samples), + m_clusters(max_clusters_per_frame), m_time_kernels(time_kernels) { if (n_streams_ <= 0) { throw std::invalid_argument( "ClusterFinderCUDA: n_streams must be > 0"); @@ -539,7 +544,11 @@ class ClusterFinderCUDA { } void push_pedestal_frame(NDView frame) { - m_pedestal.push(frame); + if (!m_pedestal.ready()) { + m_pedestal.add_init_frame(frame); + } else { + m_pedestal.push_ema(frame); + } m_pedestal_dirty = true; } diff --git a/include/aare/ClusterFinderCUDA_graph.hpp b/include/aare/ClusterFinderCUDA_graph.hpp index 98c24737..a2757a35 100644 --- a/include/aare/ClusterFinderCUDA_graph.hpp +++ b/include/aare/ClusterFinderCUDA_graph.hpp @@ -94,7 +94,7 @@ class ClusterFinderCUDAGraph { size_t m_clusters_offset; // offset of cluster array within output block COMPUTE_TYPE m_nSigma; - Pedestal m_pedestal; + FastPedestal m_pedestal; ClusterVector m_clusters; bool m_pedestal_dirty = true; bool m_graphs_dirty = true; // set when pedestal or h_output_pinned changes @@ -125,11 +125,13 @@ class ClusterFinderCUDAGraph { */ ClusterFinderCUDAGraph(Shape<2> shape_, COMPUTE_TYPE nSigma = 5.0, size_t max_clusters_per_frame = 2048, - int n_streams_ = 4) + int n_streams_ = 4, + size_t min_pedestal_samples = 1000) : m_shape(shape_), nrows(shape_[0]), ncols(shape_[1]), m_image_size(nrows * ncols), n_streams(n_streams_), m_max_clusters_per_frame(max_clusters_per_frame), m_nSigma(nSigma), - m_pedestal(shape_[0], shape_[1]), m_clusters(max_clusters_per_frame) { + m_pedestal(shape_[0], shape_[1], min_pedestal_samples), + m_clusters(max_clusters_per_frame) { if (n_streams_ <= 0) { throw std::invalid_argument( "ClusterFinderCUDAGraph: n_streams must be > 0"); @@ -268,7 +270,11 @@ class ClusterFinderCUDAGraph { } void push_pedestal_frame(NDView frame) { - m_pedestal.push(frame); + if (!m_pedestal.ready()) { + m_pedestal.add_init_frame(frame); + } else { + m_pedestal.push_ema(frame); + } m_pedestal_dirty = true; } @@ -427,8 +433,10 @@ class ClusterFinderCUDAGraph { */ void sync_pedestal_to_device() { NDArray h_mean = m_pedestal.mean(); - NDArray h_sum = m_pedestal.get_sum(); - NDArray h_sum2 = m_pedestal.get_sum2(); + // FastPedestal keeps both moments in double whatever PEDESTAL_TYPE is, + // so let the accessors name their own type here. + auto h_sum = m_pedestal.get_sum(); + auto h_sum2 = m_pedestal.get_sum2(); using DPT = device::DEVICE_PED_TYPE; const double n = static_cast(m_pedestal.n_samples()); diff --git a/include/aare/ClusterFinderFrozen.hpp b/include/aare/ClusterFinderFrozen.hpp index e9487b3a..455c77fc 100644 --- a/include/aare/ClusterFinderFrozen.hpp +++ b/include/aare/ClusterFinderFrozen.hpp @@ -4,9 +4,9 @@ #include "aare/ClusterFinder.hpp" // no_2x2_cluster guard, reused here #include "aare/ClusterVector.hpp" #include "aare/Dtype.hpp" +#include "aare/FastPedestal.hpp" #include "aare/NDArray.hpp" #include "aare/NDView.hpp" -#include "aare/Pedestal.hpp" #include "aare/defs.hpp" #include #include @@ -45,7 +45,7 @@ class ClusterFinderFrozen { PEDESTAL_TYPE m_nSigma; const PEDESTAL_TYPE c2; const PEDESTAL_TYPE c3; - Pedestal m_pedestal; + FastPedestal m_pedestal; ClusterVector m_clusters; static const uint8_t ClusterSizeX = ClusterType::cluster_size_x; @@ -54,21 +54,28 @@ class ClusterFinderFrozen { public: ClusterFinderFrozen(Shape<2> image_size, PEDESTAL_TYPE nSigma = 5.0, - size_t capacity = 1000000) + size_t capacity = 1000000, + size_t min_pedestal_samples = 1000) : m_image_size(image_size), m_nSigma(nSigma), c2(sqrt((ClusterSizeY + 1) / 2 * (ClusterSizeX + 1) / 2)), c3(sqrt(ClusterSizeX * ClusterSizeY)), - m_pedestal(image_size[0], image_size[1]), m_clusters(capacity) { + m_pedestal(image_size[0], image_size[1], min_pedestal_samples), + m_clusters(capacity) { LOG(logDEBUG) << "ClusterFinderFrozen: " << "image_size: " << image_size[0] << "x" << image_size[1] - << ", nSigma: " << nSigma << ", capacity: " << capacity; + << ", nSigma: " << nSigma << ", capacity: " << capacity + << ", min_pedestal_samples: " << min_pedestal_samples; } void set_nSigma(PEDESTAL_TYPE nSigma) { m_nSigma = nSigma; } PEDESTAL_TYPE get_nSigma() const { return m_nSigma; } void push_pedestal_frame(NDView frame) { - m_pedestal.push(frame); + if (!m_pedestal.ready()) { + m_pedestal.add_init_frame(frame); + } else { + m_pedestal.push_ema(frame); + } } NDArray pedestal() { return m_pedestal.mean(); } @@ -86,6 +93,11 @@ class ClusterFinderFrozen { } void find_clusters(NDView frame, uint64_t frame_number = 0) { + if (!m_pedestal.ready()) { + throw std::runtime_error( + "Pedestal is not ready, cannot find clusters"); + } + int dy = ClusterSizeY / 2; int dx = ClusterSizeX / 2; int has_center_pixel_x = ClusterSizeX % 2; @@ -177,7 +189,7 @@ class ClusterFinderFrozen { // only affects decisions from the NEXT frame on -> matches CUDA. Each // pixel is visited at most once per frame, so push order is irrelevant. for (const auto &p : deferred) - m_pedestal.push_fast(p.first, p.second, frame(p.first, p.second)); + m_pedestal.push_ema(p.first, p.second, frame(p.first, p.second)); } }; diff --git a/include/aare/ClusterFinderMT.hpp b/include/aare/ClusterFinderMT.hpp index 7f1f3e01..85bf0bd2 100644 --- a/include/aare/ClusterFinderMT.hpp +++ b/include/aare/ClusterFinderMT.hpp @@ -145,13 +145,12 @@ class ClusterFinderMT { * the sink */ // void collect() { - // bool empty = true; // Backoff backoff; - // while (!m_stop_requested || !empty || !m_processing_threads_stopped) - // { + // while (!m_stop_requested || !output_queues_are_empty() || + // !m_processing_threads_stopped) { // bool moved_any = false; // for (auto &queue : m_output_queues) { - // while (auto *front = queue->frontPtr()) { + // if (auto *front = queue->frontPtr(); front != nullptr) { // while (!m_sink.write(std::move(*front))) { // backoff.pause(); // } @@ -159,7 +158,6 @@ class ClusterFinderMT { // moved_any = true; // } // } - // empty = !moved_any; // if (moved_any) { // backoff.reset(); // } else { diff --git a/include/aare/FastPedestal.hpp b/include/aare/FastPedestal.hpp index 1a3579e8..0f9f1284 100644 --- a/include/aare/FastPedestal.hpp +++ b/include/aare/FastPedestal.hpp @@ -144,6 +144,43 @@ template class FastPedestal { return res; } + /** + * @brief Return a copy of the per-pixel first moment, ~ n * E[X]. + * @throws std::runtime_error if ready() is false. + * @note The two moments are stored interleaved for cache locality, so this + * de-interleaves them into a plain array. Intended for callers that need a + * contiguous moment array rather than per-pixel access, such as the CUDA + * pedestal upload. + */ + NDArray get_sum() const { + if (!ready()) { + throw std::runtime_error( + "Pedestal is not ready, cannot return sum"); + } + NDArray res({m_rows, m_cols}); + for (ssize_t i = 0; i < m_sum.size(); ++i) { + res[i] = m_sum[i].sum; + } + return res; + } + + /** + * @brief Return a copy of the per-pixel second moment, ~ n * E[X^2]. + * @throws std::runtime_error if ready() is false. + * @see get_sum() for why this is a de-interleaving copy. + */ + NDArray get_sum2() const { + if (!ready()) { + throw std::runtime_error( + "Pedestal is not ready, cannot return sum2"); + } + NDArray res({m_rows, m_cols}); + for (ssize_t i = 0; i < m_sum.size(); ++i) { + res[i] = m_sum[i].sum2; + } + return res; + } + /** * @brief Calculate the population standard deviation at (row, col). * @throws std::runtime_error if ready() is false or either index is out of diff --git a/include/aare/clusterfinder_algo.cuh b/include/aare/clusterfinder_algo.cuh index 7c5a0ae9..95ad19d4 100644 --- a/include/aare/clusterfinder_algo.cuh +++ b/include/aare/clusterfinder_algo.cuh @@ -1,8 +1,8 @@ // SPDX-License-Identifier: MPL-2.0 #pragma once #include "aare/Cluster.hpp" +#include "aare/FastPedestal.hpp" #include "aare/NDArray.hpp" -#include "aare/Pedestal.hpp" #include "aare/clusterfinder_kernel.cuh" #include "aare/utils/cuda_check.cuh" #include @@ -33,7 +33,7 @@ namespace aare::cuda { * HostImage host image of Buffers, ready to upload * Buffers allocate(size_t n_pixels) * View make_view(const Buffers &) - * HostImage prepare_pedestal(Pedestal &) [template on HostPedT] + * HostImage prepare_pedestal(FastPedestal &) [template on HostPedT] * void upload_pedestal(const HostImage &, const Buffers &, cudaStream_t) * void download_mean(const Buffers &, size_t, std::vector &, * cudaStream_t) @@ -169,10 +169,12 @@ struct CenteredRunningPedestal { * term that would not survive. */ template - static HostImage prepare_pedestal(Pedestal &host) { + static HostImage prepare_pedestal(FastPedestal &host) { NDArray h_mean = host.mean(); - NDArray h_sum = host.get_sum(); - NDArray h_sum2 = host.get_sum2(); + // FastPedestal keeps both moments in double whatever HostPedT is, so + // let the accessors name their own type here. + auto h_sum = host.get_sum(); + auto h_sum2 = host.get_sum2(); const size_t n_pixels = static_cast(h_mean.size()); const double n = static_cast(host.n_samples()); diff --git a/python/src/bind_ClusterFinderCUDA.hpp b/python/src/bind_ClusterFinderCUDA.hpp index ee8071f7..4529e61f 100644 --- a/python/src/bind_ClusterFinderCUDA.hpp +++ b/python/src/bind_ClusterFinderCUDA.hpp @@ -2,8 +2,8 @@ #pragma once #include "aare/ClusterFinderCUDA.hpp" #include "aare/ClusterVector.hpp" +#include "aare/FastPedestal.hpp" #include "aare/NDView.hpp" -#include "aare/Pedestal.hpp" #include "np_helper.hpp" #include @@ -112,9 +112,10 @@ void define_ClusterFinderCUDA(py::module &m, const std::string &typestr) { coordinates as (x, y).)"); py::class_(m, class_name.c_str()) - .def(py::init, float, size_t, int, bool>(), + .def(py::init, float, size_t, int, size_t, bool>(), py::arg("image_size"), py::arg("n_sigma") = 5.0f, py::arg("max_clusters_per_frame") = 2048, py::arg("n_streams") = 4, + py::arg("min_pedestal_samples") = 1000, py::arg("time_kernels") = false) .def_property( diff --git a/python/src/bind_ClusterFinderCUDAGraph.hpp b/python/src/bind_ClusterFinderCUDAGraph.hpp index af9dcd2a..a7af6d3b 100644 --- a/python/src/bind_ClusterFinderCUDAGraph.hpp +++ b/python/src/bind_ClusterFinderCUDAGraph.hpp @@ -2,8 +2,8 @@ #pragma once #include "aare/ClusterFinderCUDA_graph.hpp" #include "aare/ClusterVector.hpp" +#include "aare/FastPedestal.hpp" #include "aare/NDView.hpp" -#include "aare/Pedestal.hpp" #include "np_helper.hpp" #include @@ -31,9 +31,10 @@ void define_ClusterFinderCUDAGraph(py::module &m, const std::string &typestr) { py::array_t; py::class_(m, class_name.c_str()) - .def(py::init, float, size_t, int>(), py::arg("image_size"), - py::arg("n_sigma") = 5.0f, - py::arg("max_clusters_per_frame") = 2048, py::arg("n_streams") = 4) + .def(py::init, float, size_t, int, size_t>(), + py::arg("image_size"), py::arg("n_sigma") = 5.0f, + py::arg("max_clusters_per_frame") = 2048, py::arg("n_streams") = 4, + py::arg("min_pedestal_samples") = 1000) .def_property( "nSigma", &CF::get_nSigma, &CF::set_nSigma, diff --git a/python/src/bind_ClusterFinderFrozen.hpp b/python/src/bind_ClusterFinderFrozen.hpp index 1f214e4d..0efe8049 100644 --- a/python/src/bind_ClusterFinderFrozen.hpp +++ b/python/src/bind_ClusterFinderFrozen.hpp @@ -27,8 +27,10 @@ void define_ClusterFinderFrozen(py::module &m, const std::string &typestr) { py::class_>( m, class_name.c_str()) - .def(py::init, pd_type, size_t>(), py::arg("image_size"), - py::arg("n_sigma") = 5.0, py::arg("capacity") = 1'000'000) + .def(py::init, pd_type, size_t, size_t>(), + py::arg("image_size"), py::arg("n_sigma") = 5.0, + py::arg("capacity") = 1'000'000, + py::arg("min_pedestal_samples") = 1000) .def_property( "nSigma", diff --git a/python/tests/ClusterFinderCUDA.ipynb b/python/tests/ClusterFinderCUDA.ipynb index e65679b9..a23a6bce 100644 --- a/python/tests/ClusterFinderCUDA.ipynb +++ b/python/tests/ClusterFinderCUDA.ipynb @@ -88,7 +88,7 @@ "# reached, which is how you size this for a new dataset.\n", "CAP = 3000\n", "\n", - "capacity = 10_000 # initial ClusterVector capacity, CPU side\n", + "capacity = 3_000 # initial ClusterVector capacity, CPU side\n", "N_BINS = 200\n", "E_MIN, E_MAX = -2, 4000 # spectrum range, ADU\n", "\n", @@ -113,12 +113,12 @@ "outputs": [], "source": [ "cf_cpu = ClusterFinderMT(image_size, cluster_size, n_sigma=N_SIGMA,\n", - " capacity=capacity, n_threads=N_THREADS)\n", + " capacity=capacity, min_pedestal_samples=n_frames_pd, n_threads=N_THREADS)\n", "sink = ClusterCollector(cf_cpu)\n", "\n", "cf_cuda = ClusterFinderCUDA(image_size, cluster_size, n_sigma=N_SIGMA,\n", " max_clusters_per_frame=CAP, n_streams=N_STREAMS, \n", - " time_kernels=False, algorithm=\"fixed_window\")" + " n_pedestal_frames=n_frames_pd, algorithm=\"fixed_window\", time_kernels=False)" ] }, { @@ -175,7 +175,7 @@ "name": "stdout", "output_type": "stream", "text": [ - "Reading 100000 frames: 82.39s (1,214 FPS, 32.00 GB at 0.39 GB/s)\n" + "Reading 100000 frames: 5.20s (19,219 FPS, 32.00 GB at 6.15 GB/s)\n" ] } ], @@ -207,14 +207,14 @@ "name": "stderr", "output_type": "stream", "text": [ - "100%|██████████████████████████████████████████████████| 100000/100000 [00:23<00:00, 4191.14it/s]\n" + "100%|██████████████████████████████████████████████████████████████| 100000/100000 [00:23<00:00, 4310.91it/s]\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ - "CPU MT (24 threads): 23.86s 4,190 FPS 233,099,227 clusters (2331.0/frame)\n" + "CPU MT (24 threads): 23.20s 4,310 FPS 233,099,227 clusters (2331.0/frame)\n" ] } ], @@ -253,10 +253,21 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 7, "id": "cuda", "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "CUDA batched: 3.18s 31,484 FPS 233,093,554 clusters (2330.9/frame)\n", + " input pinned in 0.84s, once, outside the loop\n", + " vs CPU: 7.30x\n", + " peak clusters in one frame: 2,545 of CAP 3,000\n" + ] + } + ], "source": [ "t0 = time.perf_counter()\n", "cf_cuda.register_input_buffer(data)\n", @@ -309,10 +320,10 @@ "output_type": "stream", "text": [ " s FPS us/frame clusters\n", - "CPU MT (24 thr) 23.86 4,190 238.6 233,099,227\n", - "CUDA batched 3.34 29,926 33.4 233,093,554\n", + "CPU MT (24 thr) 23.20 4,310 232.0 233,099,227\n", + "CUDA batched 3.18 31,484 31.8 233,093,554\n", "\n", - "speedup: 7.14x\n" + "speedup: 7.30x\n" ] } ], @@ -429,7 +440,7 @@ " genuine add / drop 31,395 25,871\n", "\n", "frames touched by a genuine difference: 41,569 of 100,000\n", - "[40.6s]\n" + "[40.2s]\n" ] } ],