Implement Theil-Sen's method for fitting a line to noisy data (used in bandwidth estimation).

webrtc/modules/congestion_controller/median_slope_estimator.cc

Index: webrtc/modules/congestion_controller/median_slope_estimator.cc
diff --git a/webrtc/modules/congestion_controller/median_slope_estimator.cc b/webrtc/modules/congestion_controller/median_slope_estimator.cc
new file mode 100644
index 0000000000000000000000000000000000000000..d76ad32c1bc5797158295044c8ab02dcbaa24d88
--- /dev/null
+++ b/webrtc/modules/congestion_controller/median_slope_estimator.cc
@@ -0,0 +1,97 @@
+/*
+ *  Copyright (c) 2016 The WebRTC project authors. All Rights Reserved.
+ *
+ *  Use of this source code is governed by a BSD-style license
+ *  that can be found in the LICENSE file in the root of the source
+ *  tree. An additional intellectual property rights grant can be found
+ *  in the file PATENTS.  All contributing project authors may
+ *  be found in the AUTHORS file in the root of the source tree.
+ */
+
+#include "webrtc/modules/congestion_controller/median_slope_estimator.h"
+
+#include <algorithm>
+#include <vector>
+
+#include "webrtc/base/logging.h"
+#include "webrtc/modules/remote_bitrate_estimator/include/bwe_defines.h"
+#include "webrtc/modules/remote_bitrate_estimator/test/bwe_test_logging.h"
+
+namespace webrtc {
+
+constexpr unsigned int kDeltaCounterMax = 1000;
+
+MedianSlopeEstimator::MedianSlopeEstimator(size_t window_size,
+                                           double threshold_gain)
+    : window_size_(window_size),
+      threshold_gain_(threshold_gain),
+      num_of_deltas_(0),
+      accumulated_delay_(0),
+      delay_hist_(),
+      median_filter_(0.5),
+      trendline_(0) {}
+
+MedianSlopeEstimator::~MedianSlopeEstimator() {}
+
+void MedianSlopeEstimator::Update(double recv_delta_ms,
+                                  double send_delta_ms,
+                                  double now_ms) {
+  const double delta_ms = recv_delta_ms - send_delta_ms;
+  ++num_of_deltas_;
+  if (num_of_deltas_ > kDeltaCounterMax) {
+    num_of_deltas_ = kDeltaCounterMax;
+  }
+
+  accumulated_delay_ += delta_ms;
+  BWE_TEST_LOGGING_PLOT(1, "accumulated_delay_ms", now_ms, accumulated_delay_);
+
+  // If the window is full, remove the |window_size_| - 1 slopes that belong to
+  // the oldest point.
+  if (delay_hist_.size() == window_size_) {
+    bool success = true;
+    for (double slope : delay_hist_.front().slopes) {
+      if (!median_filter_.Erase(slope)) {
+        success = false;
+        break;
+      }
+    }
+    delay_hist_.pop_front();
+    if (!success) {
+      // We're erasing previously inserted elements from a multiset, which

[magjed_webrtc, 2016/12/08 09:31:15]

Since you don't recompute the values, I don't see why you need this. I would
just simplify the logic and do:
const bool success = median_filter_.Erase(slope);
RTC_CHECK(success);

[terelius, 2016/12/08 10:47:13]

It is basically a question of whether we prefer to crash or attempt to recover
if the unexpected happens. I haven't been able to find the relevant part of the
C++ standard, but since I don't expect it to ever happen, I don't have a strong
opinion on how it is handled. Changed as you suggested.

[magjed_webrtc, 2016/12/08 11:52:30]

I still don't understand why precision or the C++ standard is relevant as long
as you store the values and don't recompute them. The double is just bits in
memory and it won't change once calculated. Am I misunderstanding something?

+      // should always succeed. However, there are some intricacies related
+      // to extended precision of floating point registers, so we still want
+      // to verify that the operation succeeds. See comment at insertion below.
+      LOG(LS_ERROR) << "Failed to erase previously inserted slope. This is "
+                       "a bug, likely related to rounding.";
+      RTC_NOTREACHED();        // Terminate the program in debug versions.
+      median_filter_.Clear();  // Reset the filter in production versions.
+      delay_hist_.clear();
+    }
+  }
+  // Add |window_size_| - 1 new slopes.
+  for (auto& old_delay : delay_hist_) {
+    if (now_ms - old_delay.time != 0) {
+      // The C99 standard explicitly states that casts and assignments must
+      // perform the associated conversions. This means that |slope| will be
+      // a 64-bit double even if the division is computed using, e.g., 80-bit
+      // extended precision. I believe this also holds in C++ even though the
+      // C++11 standard isn't as explicit. Furthermore, there are good reasons
+      // to believe that compilers couldn't perform optimizations that break
+      // this assumption even if they wanted to.
+      double slope =
+          (accumulated_delay_ - old_delay.delay) / (now_ms - old_delay.time);
+      median_filter_.Insert(slope);
+      // We want to avoid issues with different rounding mode / precision
+      // which we might get if we recomputed the slope when we remove it.
+      old_delay.slopes.push_back(slope);
+    }
+  }
+  delay_hist_.emplace_back(now_ms, accumulated_delay_, window_size_ - 1);
+  // Recompute the median slope.
+  if (delay_hist_.size() == window_size_)
+    trendline_ = median_filter_.GetPercentileValue();
+
+  BWE_TEST_LOGGING_PLOT(1, "trendline_slope", now_ms, trendline_);
+}
+
+}  // namespace webrtc