diff --git a/src/internal_modules/roc_core/mov_quantile.h b/src/internal_modules/roc_core/mov_quantile.h index 2cdcf3f49..31620cba6 100644 --- a/src/internal_modules/roc_core/mov_quantile.h +++ b/src/internal_modules/roc_core/mov_quantile.h @@ -1,5 +1,5 @@ /* - * Copyright (c) 2023 Roc Streaming authors + * Copyright (c) 2024 Roc Streaming authors * * This Source Code Form is subject to the terms of the Mozilla Public * License, v. 2.0. If a copy of the MPL was not distributed with this @@ -19,16 +19,16 @@ namespace roc { namespace core { -//! Rolling window moving quantile -//!@remarks +//! Rolling window moving quantile. +//! @tparam T defines a sample type. +//! @remarks //! Efficiently implements moving quantile using partition heap based on approach //! described in https://aakinshin.net/posts/partitioning-heaps-quantile-estimator/. -//! It follows the quantile estimator strategy mentioned in the document -//! @tparam T defines a sample type. +//! It follows the quantile estimator strategy mentioned in the document template class MovQuantile { public: + //! Initialize MovQuantile(IArena& arena, const size_t win_len, const double quantile) - //! Initialize : win_len_(win_len) , quantile_(quantile) , old_heap_root_index_(0) @@ -40,8 +40,8 @@ template class MovQuantile { , win_filled_(false) , valid_(false) , heap_(arena) - , elem_index_heap_index_(arena) - , heap_index_elem_index_(arena) { + , elem_index_2_heap_index_(arena) + , heap_index_2_elem_index_(arena) { if (win_len == 0) { roc_panic("mov quantile: window length must be greater than 0"); } @@ -51,14 +51,14 @@ template class MovQuantile { if (!heap_.resize(win_len)) { return; } - if (!elem_index_heap_index_.resize(win_len)) { + if (!elem_index_2_heap_index_.resize(win_len)) { return; } - if (!heap_index_elem_index_.resize(win_len)) { + if (!heap_index_2_elem_index_.resize(win_len)) { return; } - double index = (quantile * static_cast(win_len - 1)); - heap_root_ = static_cast(index); + const double index = (quantile * double(win_len - 1)); + heap_root_ = size_t(index); max_heap_index_ = heap_root_; min_heap_index_ = heap_root_; valid_ = true; @@ -69,163 +69,173 @@ template class MovQuantile { return valid_; } - //! Swaps 2 heap elements along with their mapping in element to heap index and heap - //! to element index - void swap(size_t index_1, size_t index_2) { - size_t elem_index_1 = heap_index_elem_index_[index_1]; - size_t elem_index_2 = heap_index_elem_index_[index_2]; + //! Returns the moving quantile + T mov_quantile() const { + return heap_[heap_root_]; + } - T temp = heap_[index_1]; - heap_[index_1] = heap_[index_2]; - heap_[index_2] = temp; + //! Insert or swaps elements in the partition heap + //! @remarks + //! Case 1: The window is filled. The element in heap is changed whose + //! element_index%window_length is equal to arrived element. heapify_() is called post + //! that Case 2: The window in not filled. In this case we insert element in max_heap, + //! min_heap or root based on the current percentile index + void add(const T& x) { + if (elem_index_ == win_len_) { + win_filled_ = true; + } + heap_size_ = elem_index_ + 1; + elem_index_ = (elem_index_) % win_len_; + if (win_filled_) { + heap_size_ = win_len_; + min_heap_index_ = win_len_ - 1; + max_heap_index_ = 0; + const size_t heap_index = elem_index_2_heap_index_[elem_index_]; + heap_[heap_index] = x; + heapify_(heap_index); + } else { + const double index = quantile_ * double(heap_size_ - 1); + const size_t k = size_t(index); + size_t heap_index; + if (elem_index_ == 0) { + heap_index = heap_root_; + elem_index_2_heap_index_[elem_index_] = heap_index; + heap_[heap_index] = x; + heap_index_2_elem_index_[heap_index] = elem_index_; + } else { + if (old_heap_root_index_ == k) { + min_heap_index_ += 1; + heap_index = min_heap_index_; + } else { + max_heap_index_ -= 1; + heap_index = max_heap_index_; + } + elem_index_2_heap_index_[elem_index_] = heap_index; + heap_[heap_index] = x; + heap_index_2_elem_index_[heap_index] = elem_index_; + heapify_(heap_index); + old_heap_root_index_ = k; + } + } - heap_index_elem_index_[index_1] = elem_index_2; - heap_index_elem_index_[index_2] = elem_index_1; + elem_index_ = elem_index_ + 1; + } - elem_index_heap_index_[elem_index_1] = index_2; - elem_index_heap_index_[elem_index_2] = index_1; +private: + //! Maintains property of the partition heap when an element in inserted or swapped. + //! @remarks + //! The element could be inserted or changed in min_heap, max_heap or the root. + void heapify_(const size_t heap_index) { + if (heap_index < heap_root_) { + const size_t parent = heap_root_ - ((heap_root_ - heap_index - 1) / 2); + if (heap_[parent] < heap_[heap_index]) { + max_heapify_up_(heap_index); + min_heapify_down_(heap_root_); + } else { + max_heapify_down_(heap_index); + } + } else if (heap_root_ == heap_index) { + max_heapify_down_(heap_index); + min_heapify_down_(heap_root_); + } else { + const size_t parent = (heap_index - heap_root_ - 1) / 2 + heap_root_; + if (heap_[parent] > heap_[heap_index]) { + min_heapify_up_(heap_index); + max_heapify_down_(heap_root_); + } else { + min_heapify_down_(heap_index); + } + } } //! Recursively swaps parent and element in min heap partition until the parent is //! smaller or element reaches root index. - void min_heapify_up(size_t heap_index) { + void min_heapify_up_(const size_t heap_index) { if (heap_index == heap_root_) { return; } - size_t parent = (heap_index - heap_root_ - 1) / 2 + heap_root_; + const size_t parent = (heap_index - heap_root_ - 1) / 2 + heap_root_; if (heap_[parent] > heap_[heap_index]) { - swap(heap_index, parent); - min_heapify_up(parent); + swap_(heap_index, parent); + min_heapify_up_(parent); } } + //! Recursively swaps parent and element in max heap partition until the parent is //! larger or element reaches root index. //! @remarks //! The root index in max heap partition is larger than all its child index so parent //! index formulae has been adjusted accordingly - void max_heapify_up(size_t heap_index) { + void max_heapify_up_(const size_t heap_index) { // sift up if (heap_index == heap_root_) { return; } - size_t parent = heap_root_ - ((heap_root_ - heap_index - 1) / 2); + const size_t parent = heap_root_ - ((heap_root_ - heap_index - 1) / 2); if (heap_[parent] < heap_[heap_index]) { - swap(heap_index, parent); - max_heapify_up(parent); + swap_(heap_index, parent); + max_heapify_up_(parent); } } + //! Recursively swaps children and element in min heap partition until the children //! are smaller or there are no children - void min_heapify_down(size_t heap_index) { + void min_heapify_down_(const size_t heap_index) { size_t largest = heap_index; - size_t left = 2 * (heap_index - heap_root_) + 1 + heap_root_; - if (left <= min_heap_index_ && heap_[left] < heap_[largest]) + const size_t left = 2 * (heap_index - heap_root_) + 1 + heap_root_; + if (left <= min_heap_index_ && heap_[left] < heap_[largest]) { largest = left; - size_t right = 2 * (heap_index - heap_root_) + 2 + heap_root_; - if (right <= min_heap_index_ && heap_[right] < heap_[largest]) + } + const size_t right = 2 * (heap_index - heap_root_) + 2 + heap_root_; + if (right <= min_heap_index_ && heap_[right] < heap_[largest]) { largest = right; + } if (largest != heap_index) { - swap(heap_index, largest); - min_heapify_down(largest); + swap_(heap_index, largest); + min_heapify_down_(largest); } } + //! Recursively swaps children and element in max heap partition until the children //! are larger or there are no children. //! @remarks //! Similar adjustment to child index calculation like in max_heapify_up - void max_heapify_down(size_t heap_index) { + void max_heapify_down_(const size_t heap_index) { size_t largest = heap_index; - size_t left = 2 * (heap_root_ - heap_index) + 1; + const size_t left = 2 * (heap_root_ - heap_index) + 1; if (left <= (heap_root_ - max_heap_index_) - && heap_[heap_root_ - left] > heap_[largest]) + && heap_[heap_root_ - left] > heap_[largest]) { largest = heap_root_ - left; - size_t right = 2 * (heap_root_ - heap_index) + 2; + } + const size_t right = 2 * (heap_root_ - heap_index) + 2; if (right <= (heap_root_ - max_heap_index_) - && heap_[heap_root_ - right] > heap_[largest]) + && heap_[heap_root_ - right] > heap_[largest]) { largest = heap_root_ - right; - - if (largest != heap_index) { - swap(heap_index, largest); - max_heapify_down(largest); } - } - //! Maintains property of the partition heap when an element in inserted or swapped. - //! @remarks - //! The element could be inserted or changed in min_heap, max_heap or the root. - void heapify(size_t heap_index) { - if (heap_index < heap_root_) { - size_t parent = heap_root_ - ((heap_root_ - heap_index - 1) / 2); - if (heap_[parent] < heap_[heap_index]) { - max_heapify_up(heap_index); - min_heapify_down(heap_root_); - } else { - max_heapify_down(heap_index); - } - } else if (heap_root_ == heap_index) { - max_heapify_down(heap_index); - min_heapify_down(heap_root_); - } else { - size_t parent = (heap_index - heap_root_ - 1) / 2 + heap_root_; - if (heap_[parent] > heap_[heap_index]) { - min_heapify_up(heap_index); - max_heapify_down(heap_root_); - } else { - min_heapify_down(heap_index); - } + if (largest != heap_index) { + swap_(heap_index, largest); + max_heapify_down_(largest); } } - //! Insert or swaps elements in the partition heap - //! @remarks - //! Case 1: The window is filled. The element in heap is changed whose - //! element_index%window_length is equal to arrived element. heapify is called post - //! that Case 2: The window in not filled. In this case we insert element in max_heap, - //! min_heap or root based on the current percentile index - void add(const T& x) { - if (elem_index_ == win_len_) - win_filled_ = true; - heap_size_ = elem_index_ + 1; - elem_index_ = (elem_index_) % win_len_; - if (win_filled_) { - heap_size_ = win_len_; - min_heap_index_ = win_len_ - 1; - max_heap_index_ = 0; - size_t heap_index = elem_index_heap_index_[elem_index_]; - heap_[heap_index] = x; - heapify(heap_index); - } else { - double index = quantile_ * static_cast(heap_size_ - 1); - size_t k = static_cast(index); - size_t heap_index; - if (elem_index_ == 0) { - heap_index = heap_root_; - elem_index_heap_index_[elem_index_] = heap_index; - heap_[heap_index] = x; - heap_index_elem_index_[heap_index] = elem_index_; - } else { - if (old_heap_root_index_ == k) { - min_heap_index_ += 1; - heap_index = min_heap_index_; - } else { - max_heap_index_ -= 1; - heap_index = max_heap_index_; - } - elem_index_heap_index_[elem_index_] = heap_index; - heap_[heap_index] = x; - heap_index_elem_index_[heap_index] = elem_index_; - heapify(heap_index); - old_heap_root_index_ = k; - } - } + //! Swaps 2 heap elements along with their mapping in element to heap index and heap + //! to element index + void swap_(const size_t index_1, const size_t index_2) { + const size_t elem_index_1 = heap_index_2_elem_index_[index_1]; + const size_t elem_index_2 = heap_index_2_elem_index_[index_2]; - elem_index_ = elem_index_ + 1; - } - //! Returns the moving quantile - T sliding_quantile() { - return heap_[heap_root_]; + const T temp = heap_[index_1]; + heap_[index_1] = heap_[index_2]; + heap_[index_2] = temp; + + heap_index_2_elem_index_[index_1] = elem_index_2; + heap_index_2_elem_index_[index_2] = elem_index_1; + + elem_index_2_heap_index_[elem_index_1] = index_2; + elem_index_2_heap_index_[elem_index_2] = index_1; } //! Length of the sliding window @@ -235,30 +245,30 @@ template class MovQuantile { //! Used to check the window filling logic size_t old_heap_root_index_; - //! Index which seperates max and min heap and also act as their root + //! Index which separates max and min heap and also act as their root size_t heap_root_; //! Maintains current heap size size_t heap_size_; - //! Maintians the index to which max_heap extends + //! Maintains the index to which max_heap extends size_t max_heap_index_; - //! Maintians the index to which min_heap extends + //! Maintains the index to which min_heap extends size_t min_heap_index_; - //! Maintians current element index + //! Maintains current element index size_t elem_index_; //! Window filled check bool win_filled_; - //! Maintians initialization success + //! Maintains initialization success bool valid_; - //! Maintians the parititon heap + //! Maintains the partition heap Array heap_; - //! Maintians the element index to heap index mapping - Array elem_index_heap_index_; + //! Maintains the element index to heap index mapping + Array elem_index_2_heap_index_; //! Maintains the heap index to element index mapping - Array heap_index_elem_index_; + Array heap_index_2_elem_index_; }; } // namespace core diff --git a/src/internal_modules/roc_core/mov_stats.h b/src/internal_modules/roc_core/mov_stats.h index b2dc9ec97..37b40c906 100644 --- a/src/internal_modules/roc_core/mov_stats.h +++ b/src/internal_modules/roc_core/mov_stats.h @@ -21,11 +21,10 @@ namespace roc { namespace core { //! Rolling window moving average and variance. -//! @remarks -//! Efficiently implements moving average and variance based on approach -//! described in https://www.dsprelated.com/showthread/comp.dsp/97276-1.php -//! //! @tparam T defines a sample type. +//! @remarks +//! Efficiently implements moving average and variance based on approach +//! described in https://www.dsprelated.com/showthread/comp.dsp/97276-1.php template class MovStats { public: //! Initialize. diff --git a/src/tests/roc_core/test_mov_quantile.cpp b/src/tests/roc_core/test_mov_quantile.cpp index fd1cad78d..552b9d0db 100644 --- a/src/tests/roc_core/test_mov_quantile.cpp +++ b/src/tests/roc_core/test_mov_quantile.cpp @@ -1,5 +1,5 @@ /* - * Copyright (c) 2023 Roc Streaming authors + * Copyright (c) 2024 Roc Streaming authors * * This Source Code Form is subject to the terms of the Mozilla Public * License, v. 2.0. If a copy of the MPL was not distributed with this @@ -14,11 +14,11 @@ namespace roc { namespace core { -TEST_GROUP(movquantile) { +TEST_GROUP(mov_quantile) { HeapArena arena; }; -TEST(movquantile, testing_minimum) { +TEST(mov_quantile, testing_minimum) { const size_t n = 9; MovQuantile quant(arena, n, 0.0); CHECK(quant.is_valid()); @@ -29,16 +29,16 @@ TEST(movquantile, testing_minimum) { quant.add(18); quant.add(15); quant.add(25); - LONGS_EQUAL((int64_t)11, quant.sliding_quantile()); // test window incomplete + LONGS_EQUAL((int64_t)11, quant.mov_quantile()); // test window incomplete quant.add(32); quant.add(14); quant.add(19); quant.add(16); quant.add(35); - LONGS_EQUAL((int64_t)12, quant.sliding_quantile()); // test window complete + LONGS_EQUAL((int64_t)12, quant.mov_quantile()); // test window complete } -TEST(movquantile, testing_lower_side) { +TEST(mov_quantile, testing_lower_side) { const size_t n = 12; MovQuantile quant(arena, n, 0.34); CHECK(quant.is_valid()); @@ -49,7 +49,7 @@ TEST(movquantile, testing_lower_side) { quant.add(18); quant.add(6); quant.add(24); - LONGS_EQUAL((int64_t)12, quant.sliding_quantile()); // test window incomplete + LONGS_EQUAL((int64_t)12, quant.mov_quantile()); // test window incomplete quant.add(22); quant.add(35); quant.add(42); @@ -60,11 +60,11 @@ TEST(movquantile, testing_lower_side) { quant.add(45); quant.add(49); quant.add(37); - int64_t x1 = quant.sliding_quantile(); // test complete window insertion + int64_t x1 = quant.mov_quantile(); // test complete window insertion LONGS_EQUAL((int64_t)24, x1); } -TEST(movquantile, testing_median) { +TEST(mov_quantile, testing_median) { const size_t n = 10; MovQuantile quant(arena, n, 0.50); CHECK(quant.is_valid()); @@ -75,7 +75,7 @@ TEST(movquantile, testing_median) { quant.add(25); quant.add(6); quant.add(37); - LONGS_EQUAL((int64_t)25, quant.sliding_quantile()); // test window incomplete + LONGS_EQUAL((int64_t)25, quant.mov_quantile()); // test window incomplete quant.add(23); quant.add(48); quant.add(100); @@ -86,10 +86,10 @@ TEST(movquantile, testing_median) { quant.add(72); quant.add(83); quant.add(37); - LONGS_EQUAL((int64_t)57, quant.sliding_quantile()); // test complete window + LONGS_EQUAL((int64_t)57, quant.mov_quantile()); // test complete window } -TEST(movquantile, testing_upper_side) { +TEST(mov_quantile, testing_upper_side) { const size_t n = 11; MovQuantile quant(arena, n, 0.78); CHECK(quant.is_valid()); @@ -101,16 +101,16 @@ TEST(movquantile, testing_upper_side) { quant.add(52); quant.add(14); quant.add(46); - LONGS_EQUAL((int64_t)39, quant.sliding_quantile()); // test incomplete window + LONGS_EQUAL((int64_t)39, quant.mov_quantile()); // test incomplete window quant.add(14); quant.add(14); quant.add(100); quant.add(32); quant.add(83); - LONGS_EQUAL((int64_t)46, quant.sliding_quantile()); // test complete window + LONGS_EQUAL((int64_t)46, quant.mov_quantile()); // test complete window } -TEST(movquantile, test_maximum) { +TEST(mov_quantile, test_maximum) { const size_t n = 7; MovQuantile quant(arena, n, 1); CHECK(quant.is_valid()); @@ -119,7 +119,7 @@ TEST(movquantile, test_maximum) { quant.add(38); quant.add(72); quant.add(63); - LONGS_EQUAL((int64_t)72, quant.sliding_quantile()); // test incomplete window + LONGS_EQUAL((int64_t)72, quant.mov_quantile()); // test incomplete window quant.add(35); quant.add(76); quant.add(42); @@ -129,7 +129,7 @@ TEST(movquantile, test_maximum) { quant.add(102); quant.add(56); quant.add(20); - LONGS_EQUAL((int64_t)102, quant.sliding_quantile()); // test complete window + LONGS_EQUAL((int64_t)102, quant.mov_quantile()); // test complete window } } // namespace core diff --git a/src/tests/roc_core/test_mov_stats.cpp b/src/tests/roc_core/test_mov_stats.cpp index 4b74c6e9f..642ecaa7b 100644 --- a/src/tests/roc_core/test_mov_stats.cpp +++ b/src/tests/roc_core/test_mov_stats.cpp @@ -42,11 +42,11 @@ long Object::n_objects = 0; } // namespace -TEST_GROUP(movstats) { +TEST_GROUP(mov_stats) { HeapArena arena; }; -TEST(movstats, single_pass) { +TEST(mov_stats, single_pass) { const size_t n = 10; int64_t x[n]; MovStats stats(arena, n); @@ -69,7 +69,7 @@ TEST(movstats, single_pass) { } } -TEST(movstats, one_n_half_pass) { +TEST(mov_stats, one_n_half_pass) { const size_t n = 10; MovStats stats(arena, n); for (size_t i = 0; i < (n * 10 + n / 2); i++) { @@ -93,7 +93,7 @@ TEST(movstats, one_n_half_pass) { LONGS_EQUAL(target_var, stats.mov_var()); } -TEST(movstats, one_n_half_extend) { +TEST(mov_stats, one_n_half_extend) { const size_t n = 10; MovStats stats(arena, n); const int64_t target_avg = n;