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

webrtc/modules/congestion_controller/theil_sen_estimator.cc

Index: webrtc/modules/congestion_controller/theil_sen_estimator.cc
diff --git a/webrtc/modules/congestion_controller/theil_sen_estimator.cc b/webrtc/modules/congestion_controller/theil_sen_estimator.cc
new file mode 100644
index 0000000000000000000000000000000000000000..ba2858c8481a159285a6652f7d4dd4896ef6696e
--- /dev/null
+++ b/webrtc/modules/congestion_controller/theil_sen_estimator.cc
@@ -0,0 +1,78 @@
+/*
+ *  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/theil_sen_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 {
+
+enum { kDeltaCounterMax = 1000 };

[brandtr, 2016/11/28 16:27:19]

Can you use a constexpr size_t instead?

[terelius, 2016/12/02 16:45:52]

Done, except using unsigned int.

+
+TheilSenEstimator::TheilSenEstimator(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) {}
+
+TheilSenEstimator::~TheilSenEstimator() {}
+
+void TheilSenEstimator::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_);
+
+  // First remove the oldest N-1 slopes if the window is full.

[brandtr, 2016/11/28 16:27:19]

This could be a computer science thing, but it's not immediately obvious to me
that |N| corresponds to the number of points in the data structure. Would be
nice to point this out, for SP people like me ;)

Also, consider rephrasing comment as: "remove the N-1 slopes that belong to the
oldest point". The idea being that just because a slope partially belongs to an
old point, doesn't mean that the slope itself is old. E.g., the slope between
the oldest point and the newest point is just as new as the slope between the
second newest point and the newest point.

[terelius, 2016/12/02 16:45:52]

Done.

+  if (delay_hist_.size() == window_size_) {
+    std::pair<double, double> old_delay = delay_hist_.front();

[stefan-webrtc, 2016/11/28 13:01:58]

I think auto would be fine here if you'd prefer that.

[terelius, 2016/12/02 16:45:52]

This part of the code has changed, so the comment no longer applies.

+    delay_hist_.pop_front();
+    for (const auto& new_delay : delay_hist_) {
+      if (new_delay.first - old_delay.first != 0) {
+        double slope = (new_delay.second - old_delay.second) /
+                       (new_delay.first - old_delay.first);

[brandtr, 2016/11/28 16:27:19]

I see the problem here that you were discussing before.

One possible solution, which would not involve changing the percentile filter
interface, would be to quantize the slopes before inserting them into the
filter. As long as the quantization error is larger than the numerical error in
calculating the division, you should always find the right slope in the filter.
And as long as the quantization error is smaller than the natural variation of
delay deltas, the quantization should not impact the accuracy of the median
calculation in any noteworthy way. After quantization, you could represent the
slopes using integers.

For millisecond samples, I guess quantizing to the nearest 1e-3 ~ 1e-6 would be
enough. I don't know how large your numerical errors are though?

[ I believe that the floating points are not evenly spaced along the number
line. I don't know if that fact impacts my arguments here, as I have been
thinking in terms of real numbers :) ]

[terelius, 2016/12/02 16:45:52]

A problem with the quantization approach is that I don't a priori know the
interval/distribution of the data, nor the maximum numerical error I can get
from fp registers larger than 64 bits.

[brandtr, 2016/12/05 12:34:44]

Right. This comment, as well as the one on line 63, assumes that we will never
encounter infinitely dynamic channels; or in other words, that we have a
"bandwidth" (in the EE sense) that is bounded. But as you explained offline,
this assumption doesn't hold due to how the slopes are defined.

+        if (!median_filter_.Erase(slope))

[stefan-webrtc, 2016/11/28 13:01:58]

{}

[terelius, 2016/12/02 16:45:52]

Done.

+          LOG(LS_ERROR) << "Failed to erase previously inserted slope. This is "
+                           "a bug, likely related to rounding.";

[stefan-webrtc, 2016/11/28 13:01:58]

Should this happen? Otherwise a DCHECK might make sense?

[terelius, 2016/12/02 16:45:52]

Kept the error message but added an RTC_NOTREACHED();

+      }
+    }
+  }
+  // Add N-1 new slopes.
+  for (const auto& old_delay : delay_hist_) {
+    if (now_ms - old_delay.first != 0) {

[brandtr, 2016/11/28 16:27:19]

I guess that if (now_ms - old_delay.first) is close to 0, you will get some
unrealistic slope calculations. Would this happen on mobile networks with large
buffers and channel-dependent scheduling? If not, it could be an idea to
saturate the absolute value of the slopes at some high limit.

[terelius, 2016/12/02 16:45:52]

The filter should be robust anyway, so I don't think this warrants special
treatment. Also note that clipping some values wont affect the median unless the
median itself is outside your limit.

+      double slope =
+          (accumulated_delay_ - old_delay.second) / (now_ms - old_delay.first);
+      median_filter_.Insert(slope);
+    }
+  }
+  delay_hist_.push_back(std::make_pair(now_ms, accumulated_delay_));
+  // Recompute the median slope.
+  if (delay_hist_.size() == window_size_) {

[stefan-webrtc, 2016/11/28 13:01:58]

Feel free to remove {}

[terelius, 2016/12/02 16:45:52]

Done.

+    trendline_ = median_filter_.GetPercentileValue();
+  }
+
+  BWE_TEST_LOGGING_PLOT(1, "trendline_slope", now_ms, trendline_);
+}
+
+}  // namespace webrtc