Improve CMedianFilter algorithm

src/Makefile.bench.include

@@ -45,6 +45,7 @@ bench_bench_bitcoin_SOURCES = \
   bench/rpc_blockchain.cpp \
   bench/rpc_mempool.cpp \
   bench/strencodings.cpp \
+  bench/timedata.cpp \
   bench/util_time.cpp \
   bench/verify_script.cpp
 

src/bench/timedata.cpp

@@ -0,0 +1,72 @@
+// Copyright (c) 2020-2022 The Bitcoin Core developers
+// Distributed under the MIT software license, see the accompanying
+// file COPYING or http://www.opensource.org/licenses/mit-license.php.
+
+#include <timedata.h>
+
+#include <bench/bench.h>
+
+// 256 element max
+// 1000 iteration
+static void ComputeMedian256_1000(benchmark::Bench& bench)
+{
+    //Precompute random numbers
+    std::vector<int> rand_num;
+    rand_num.reserve(1000);
+    for (size_t i = 0; i < 1000; ++i) {
+        rand_num.emplace_back(rand());
+    }
+
+    bench.run([&rand_num] {
+        CMedianFilter<int, 256> median_filter(15);
+
+        for (size_t i = 0; i < 1000; ++i) {
+            median_filter.input(rand_num[i]);
+        }
+    });
+}
+
+// 20 element max
+// 1000 iteration
+static void ComputeMedian20_1000(benchmark::Bench& bench)
+{
+    //Precompute random numbers
+    std::vector<int> rand_num;
+    rand_num.reserve(1000);
+    for (size_t i = 0; i < 1000; ++i) {
+        rand_num.emplace_back(rand());
+    }
+
+    bench.run([&rand_num] {
+        CMedianFilter<int, 20> median_filter(15);
+
+        for (size_t i = 0; i < 1000; ++i) {
+            median_filter.input(rand_num[i]);
+        }
+    });
+}
+
+// 10 element max
+// 1000 iteration
+static void ComputeMedian10_1000(benchmark::Bench& bench)
+{
+    //Precompute random numbers
+    std::vector<int> rand_num;
+    rand_num.reserve(1000);
+    for (size_t i = 0; i < 1000; ++i) {
+        rand_num.emplace_back(rand());
+    }
+
+    bench.run([&rand_num] {
+        CMedianFilter<int, 10> median_filter(15);
+
+        for (size_t i = 0; i < 1000; ++i) {
+            median_filter.input(rand_num[i]);
+        }
+    });
+}
+
+BENCHMARK(ComputeMedian256_1000);
+BENCHMARK(ComputeMedian20_1000);
+BENCHMARK(ComputeMedian10_1000);
+

src/test/fuzz/timedata.cpp

@@ -14,12 +14,11 @@
 FUZZ_TARGET(timedata)
 {
     FuzzedDataProvider fuzzed_data_provider(buffer.data(), buffer.size());
-    const unsigned int max_size = fuzzed_data_provider.ConsumeIntegralInRange<unsigned int>(0, 1000);
-    // A max_size of 0 implies no limit, so cap the max number of insertions to avoid timeouts
+    constexpr unsigned int max_size = 200;
     auto max_to_insert = fuzzed_data_provider.ConsumeIntegralInRange<int>(0, 4000);
     // Divide by 2 to avoid signed integer overflow in .median()
     const int64_t initial_value = fuzzed_data_provider.ConsumeIntegral<int64_t>() / 2;
-    CMedianFilter<int64_t> median_filter{max_size, initial_value};
+    CMedianFilter<int64_t, max_size> median_filter(initial_value);
     while (fuzzed_data_provider.remaining_bytes() > 0 && --max_to_insert >= 0) {
         (void)median_filter.median();
         assert(median_filter.size() > 0);

src/test/timedata_tests.cpp

@@ -20,7 +20,7 @@ BOOST_FIXTURE_TEST_SUITE(timedata_tests, BasicTestingSetup)
 
 BOOST_AUTO_TEST_CASE(util_MedianFilter)
 {
-    CMedianFilter<int> filter(5, 15);
+    CMedianFilter<int, 5> filter(15);
 
     BOOST_CHECK_EQUAL(filter.median(), 15);
 
@@ -43,6 +43,87 @@ BOOST_AUTO_TEST_CASE(util_MedianFilter)
     BOOST_CHECK_EQUAL(filter.median(), 7);
 }
 
+BOOST_AUTO_TEST_CASE(util_MedianFilter2)
+{
+    CMedianFilter<int, 5> filter(0);
+
+    filter.input(2); // [0 2]
+    filter.input(5); // [0 2 5]
+    filter.input(7); // [0 2 5 7]
+    filter.input(8); // [0 2 5 7 8]
+
+    filter.input(6); // [2 5 6 7 8]
+    BOOST_CHECK_EQUAL(filter.median(), 6);
+}
+
+BOOST_AUTO_TEST_CASE(util_MedianFilter3)
+{
+    CMedianFilter<int, 5> filter(7);
+
+    filter.input(2); // [2 7]
+    filter.input(5); // [2 5 7]
+    filter.input(0); // [0 2 5 7]
+    filter.input(8); // [0 2 5 7 8]
+
+    filter.input(1); // [0 1 2 5 8]
+    BOOST_CHECK_EQUAL(filter.median(), 2);
+}
+
+BOOST_AUTO_TEST_CASE(util_MedianFilter4)
+{
+    CMedianFilter<int, 2> filter(0);
+
+    filter.input(2); // [0 2]
+    BOOST_CHECK_EQUAL(filter.median(), 1);
+
+    filter.input(3); // [2 3]
+    BOOST_CHECK_EQUAL(filter.median(), 2);
+
+    filter.input(7); // [3 7]
+    BOOST_CHECK_EQUAL(filter.median(), 5);
+
+    filter.input(7); // [7 7]
+    BOOST_CHECK_EQUAL(filter.median(), 7);
+
+    filter.input(0); // [0 7]
+    BOOST_CHECK_EQUAL(filter.median(), 3);
+}
+
+BOOST_AUTO_TEST_CASE(util_MedianFilter_rand)
+{
+    // Use a ring buffer to maintain the last 20 elements and sort them when
+    // it is time to compute the median
+    //
+    // CMedianFilter "maintains" a sorted array which should be equal to the one above
+
+    constexpr size_t ring_buffer_size = 20;
+    CMedianFilter<int, ring_buffer_size> filter(0);
+
+    std::array<int, ring_buffer_size> values;
+    size_t idx = 0;
+
+    for (size_t i = 0; i < 1'000'000; ++i) {
+        const int r = rand();
+        filter.input(r);
+
+        values[idx++] = r;
+        if (idx >= ring_buffer_size) idx = 0;
+
+        // Skip check since they won't be equal until they reach 20 elements each
+        if (i < ring_buffer_size) continue;
+
+        std::vector<int> sorted_values(values.begin(), values.end());
+        std::sort(sorted_values.begin(), sorted_values.end());
+
+        // Check that both sorted arrays are the same
+        const std::vector<int> filter_sorted_values = filter.sorted();
+        BOOST_CHECK_EQUAL(filter_sorted_values.size(), sorted_values.size());
+        for (size_t j = 0; j < sorted_values.size(); ++j) {
+            BOOST_CHECK_EQUAL(filter_sorted_values[j], sorted_values[j]);
+        }
+    }
+}
+
 static void MultiAddTimeData(int n, int64_t offset)
 {
     static int cnt = 0;

src/timedata.cpp

@@ -40,7 +40,7 @@ int64_t GetAdjustedTime()
 #define BITCOIN_TIMEDATA_MAX_SAMPLES 200
 
 static std::set<CNetAddr> g_sources;
-static CMedianFilter<int64_t> g_time_offsets{BITCOIN_TIMEDATA_MAX_SAMPLES, 0};
+static CMedianFilter<int64_t, BITCOIN_TIMEDATA_MAX_SAMPLES> g_time_offsets(0);
 static bool g_warning_emitted;
 
 void AddTimeData(const CNetAddr& ip, int64_t nOffsetSample)
@@ -75,7 +75,7 @@ void AddTimeData(const CNetAddr& ip, int64_t nOffsetSample)
     //
     if (g_time_offsets.size() >= 5 && g_time_offsets.size() % 2 == 1) {
         int64_t nMedian = g_time_offsets.median();
-        std::vector<int64_t> vSorted = g_time_offsets.sorted();
+        const std::vector<int64_t>& vSorted = g_time_offsets.sorted();
         // Only let other nodes change our time by so much
         int64_t max_adjustment = std::max<int64_t>(0, gArgs.GetIntArg("-maxtimeadjustment", DEFAULT_MAX_TIME_ADJUSTMENT));
         if (nMedian >= -max_adjustment && nMedian <= max_adjustment) {
@@ -115,6 +115,6 @@ void TestOnlyResetTimeData()
     LOCK(g_timeoffset_mutex);
     nTimeOffset = 0;
     g_sources.clear();
-    g_time_offsets = CMedianFilter<int64_t>{BITCOIN_TIMEDATA_MAX_SAMPLES, 0};
+    g_time_offsets = CMedianFilter<int64_t, BITCOIN_TIMEDATA_MAX_SAMPLES>(0);
     g_warning_emitted = false;
 }

src/timedata.h

@@ -6,7 +6,9 @@
 #define BITCOIN_TIMEDATA_H
 
 #include <algorithm>
+#include <array>
 #include <assert.h>
+#include <cstring>
 #include <stdint.h>
 #include <vector>
 
@@ -18,53 +20,75 @@ class CNetAddr;
  * Median filter over a stream of values.
  * Returns the median of the last N numbers
  */
-template <typename T>
+template <typename T, uint32_t N>
 class CMedianFilter
 {
 private:
-    std::vector<T> vValues;
     std::vector<T> vSorted;
-    unsigned int nSize;
+    std::array<T, N> vValues;
+    size_t vIdx;
 
 public:
-    CMedianFilter(unsigned int _size, T initial_value) : nSize(_size)
+    CMedianFilter(T initial_value) : vIdx(0)
     {
-        vValues.reserve(_size);
-        vValues.push_back(initial_value);
-        vSorted = vValues;
+        static_assert(std::is_trivially_copyable<T>::value);
+        static_assert(N > 1);
+        static_assert(N <= 256,
+                      "Implement algorithm that does not require maintaining a sorted array");
+
+        vValues[vIdx++] = initial_value;
+        vSorted.reserve(N);
+        vSorted.emplace_back(initial_value);
     }
 
     void input(T value)
     {
-        if (vValues.size() == nSize) {
-            vValues.erase(vValues.begin());
+        T tmp = vValues[vIdx];
+        vValues[vIdx++] = value;
+        if (vIdx >= N) vIdx = 0;
+
+        if (vSorted.size() < N) {
+            vSorted.insert(std::upper_bound(vSorted.begin(), vSorted.end(), value), value);
+            return;
         }
-        vValues.push_back(value);
 
-        vSorted.resize(vValues.size());
-        std::copy(vValues.begin(), vValues.end(), vSorted.begin());
-        std::sort(vSorted.begin(), vSorted.end());
+        // Find current element to remove and where new element should be inserted
+        // [0, 2, 2, 4, 5 ,7, 8, 9, 11, 14, 15, 16]
+        //        ^              ^
+        //        e              i
+        // Shift entire block of memory either left or right and insert new value
+
+        const auto erase_it = std::find(vSorted.begin(), vSorted.end(), tmp);
+        auto insert_it = std::upper_bound(vSorted.begin(), vSorted.end(), value);
+        assert(erase_it != vSorted.end());
+
+        if (erase_it == insert_it) {
+            *insert_it = value;
+        } else if (erase_it < insert_it) {
+            // Left rotate ensures |insert_it > vSorted.begin()|
+            std::rotate(erase_it, erase_it + 1, insert_it);
+            --insert_it;
+            *insert_it = value;
+        } else {
+            // Right rotate
+            std::rotate(insert_it, erase_it, erase_it + 1);
+            *insert_it = value;
+        }
     }
 
     T median() const
     {
-        int vSortedSize = vSorted.size();
-        assert(vSortedSize > 0);
-        if (vSortedSize & 1) // Odd number of elements
-        {
-            return vSorted[vSortedSize / 2];
-        } else // Even number of elements
-        {
-            return (vSorted[vSortedSize / 2 - 1] + vSorted[vSortedSize / 2]) / 2;
-        }
+        return vSorted.size() & 1 ?
+                   vSorted[vSorted.size() / 2] :
+                   (vSorted[vSorted.size() / 2 - 1] + vSorted[vSorted.size() / 2]) / 2;
     }
 
-    int size() const
+    size_t size() const
     {
-        return vValues.size();
+        return vSorted.size();
     }
 
-    std::vector<T> sorted() const
+    const std::vector<T>& sorted()
     {
         return vSorted;
     }