1 files changed, 206 insertions, 0 deletions
diff --git a/third_party/libwebrtc/rtc_tools/frame_analyzer/linear_least_squares.cc b/third_party/libwebrtc/rtc_tools/frame_analyzer/linear_least_squares.cc
new file mode 100644
index 0000000000..93a6f90e69
--- /dev/null
+++ b/third_party/libwebrtc/rtc_tools/frame_analyzer/linear_least_squares.cc
@@ -0,0 +1,206 @@
+/*
+ *  Copyright (c) 2018 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 "rtc_tools/frame_analyzer/linear_least_squares.h"
+
+#include <math.h>
+
+#include <cstdint>
+#include <cstdlib>
+#include <functional>
+#include <numeric>
+#include <type_traits>
+#include <utility>
+
+#include "rtc_base/checks.h"
+#include "rtc_base/logging.h"
+
+namespace webrtc {
+namespace test {
+
+template <class T>
+using Matrix = std::valarray<std::valarray<T>>;
+
+namespace {
+
+template <typename R, typename T>
+R DotProduct(const std::valarray<T>& a, const std::valarray<T>& b) {
+  RTC_CHECK_EQ(a.size(), b.size());
+  return std::inner_product(std::begin(a), std::end(a), std::begin(b), R(0));
+}
+
+// Calculates a^T * b.
+template <typename R, typename T>
+Matrix<R> MatrixMultiply(const Matrix<T>& a, const Matrix<T>& b) {
+  Matrix<R> result(std::valarray<R>(a.size()), b.size());
+  for (size_t i = 0; i < a.size(); ++i) {
+    for (size_t j = 0; j < b.size(); ++j)
+      result[j][i] = DotProduct<R>(a[i], b[j]);
+  }
+
+  return result;
+}
+
+template <typename T>
+Matrix<T> Transpose(const Matrix<T>& matrix) {
+  if (matrix.size() == 0)
+    return Matrix<T>();
+  const size_t rows = matrix.size();
+  const size_t columns = matrix[0].size();
+  Matrix<T> result(std::valarray<T>(rows), columns);
+
+  for (size_t i = 0; i < rows; ++i) {
+    for (size_t j = 0; j < columns; ++j)
+      result[j][i] = matrix[i][j];
+  }
+
+  return result;
+}
+
+// Convert valarray from type T to type R.
+template <typename R, typename T>
+std::valarray<R> ConvertTo(const std::valarray<T>& v) {
+  std::valarray<R> result(v.size());
+  for (size_t i = 0; i < v.size(); ++i)
+    result[i] = static_cast<R>(v[i]);
+  return result;
+}
+
+// Convert valarray Matrix from type T to type R.
+template <typename R, typename T>
+Matrix<R> ConvertTo(const Matrix<T>& mat) {
+  Matrix<R> result(mat.size());
+  for (size_t i = 0; i < mat.size(); ++i)
+    result[i] = ConvertTo<R>(mat[i]);
+  return result;
+}
+
+// Convert from valarray Matrix back to the more conventional std::vector.
+template <typename T>
+std::vector<std::vector<T>> ToVectorMatrix(const Matrix<T>& m) {
+  std::vector<std::vector<T>> result;
+  for (const std::valarray<T>& v : m)
+    result.emplace_back(std::begin(v), std::end(v));
+  return result;
+}
+
+// Create a valarray Matrix from a conventional std::vector.
+template <typename T>
+Matrix<T> FromVectorMatrix(const std::vector<std::vector<T>>& mat) {
+  Matrix<T> result(mat.size());
+  for (size_t i = 0; i < mat.size(); ++i)
+    result[i] = std::valarray<T>(mat[i].data(), mat[i].size());
+  return result;
+}
+
+// Returns `matrix_to_invert`^-1 * `right_hand_matrix`. `matrix_to_invert` must
+// have square size.
+Matrix<double> GaussianElimination(Matrix<double> matrix_to_invert,
+                                   Matrix<double> right_hand_matrix) {
+  // `n` is the width/height of `matrix_to_invert`.
+  const size_t n = matrix_to_invert.size();
+  // Make sure `matrix_to_invert` has square size.
+  for (const std::valarray<double>& column : matrix_to_invert)
+    RTC_CHECK_EQ(n, column.size());
+  // Make sure `right_hand_matrix` has correct size.
+  for (const std::valarray<double>& column : right_hand_matrix)
+    RTC_CHECK_EQ(n, column.size());
+
+  // Transpose the matrices before and after so that we can perform Gaussian
+  // elimination on the columns instead of the rows, since that is easier with
+  // our representation.
+  matrix_to_invert = Transpose(matrix_to_invert);
+  right_hand_matrix = Transpose(right_hand_matrix);
+
+  // Loop over the diagonal of `matrix_to_invert` and perform column reduction.
+  // Column reduction is a sequence of elementary column operations that is
+  // performed on both `matrix_to_invert` and `right_hand_matrix` until
+  // `matrix_to_invert` has been transformed to the identity matrix.
+  for (size_t diagonal_index = 0; diagonal_index < n; ++diagonal_index) {
+    // Make sure the diagonal element has the highest absolute value by
+    // swapping columns if necessary.
+    for (size_t column = diagonal_index + 1; column < n; ++column) {
+      if (std::abs(matrix_to_invert[column][diagonal_index]) >
+          std::abs(matrix_to_invert[diagonal_index][diagonal_index])) {
+        std::swap(matrix_to_invert[column], matrix_to_invert[diagonal_index]);
+        std::swap(right_hand_matrix[column], right_hand_matrix[diagonal_index]);
+      }
+    }
+
+    // Reduce the diagonal element to be 1, by dividing the column with that
+    // value. If the diagonal element is 0, it means the system of equations has
+    // many solutions, and in that case we will return an arbitrary solution.
+    if (matrix_to_invert[diagonal_index][diagonal_index] == 0.0) {
+      RTC_LOG(LS_WARNING) << "Matrix is not invertible, ignoring.";
+      continue;
+    }
+    const double diagonal_element =
+        matrix_to_invert[diagonal_index][diagonal_index];
+    matrix_to_invert[diagonal_index] /= diagonal_element;
+    right_hand_matrix[diagonal_index] /= diagonal_element;
+
+    // Eliminate the other entries in row `diagonal_index` by making them zero.
+    for (size_t column = 0; column < n; ++column) {
+      if (column == diagonal_index)
+        continue;
+      const double row_element = matrix_to_invert[column][diagonal_index];
+      matrix_to_invert[column] -=
+          row_element * matrix_to_invert[diagonal_index];
+      right_hand_matrix[column] -=
+          row_element * right_hand_matrix[diagonal_index];
+    }
+  }
+
+  // Transpose the result before returning it, explained in comment above.
+  return Transpose(right_hand_matrix);
+}
+
+}  // namespace
+
+IncrementalLinearLeastSquares::IncrementalLinearLeastSquares() = default;
+IncrementalLinearLeastSquares::~IncrementalLinearLeastSquares() = default;
+
+void IncrementalLinearLeastSquares::AddObservations(
+    const std::vector<std::vector<uint8_t>>& x,
+    const std::vector<std::vector<uint8_t>>& y) {
+  if (x.empty() || y.empty())
+    return;
+  // Make sure all columns are the same size.
+  const size_t n = x[0].size();
+  for (const std::vector<uint8_t>& column : x)
+    RTC_CHECK_EQ(n, column.size());
+  for (const std::vector<uint8_t>& column : y)
+    RTC_CHECK_EQ(n, column.size());
+
+  // We will multiply the uint8_t values together, so we need to expand to a
+  // type that can safely store those values, i.e. uint16_t.
+  const Matrix<uint16_t> unpacked_x = ConvertTo<uint16_t>(FromVectorMatrix(x));
+  const Matrix<uint16_t> unpacked_y = ConvertTo<uint16_t>(FromVectorMatrix(y));
+
+  const Matrix<uint64_t> xx = MatrixMultiply<uint64_t>(unpacked_x, unpacked_x);
+  const Matrix<uint64_t> xy = MatrixMultiply<uint64_t>(unpacked_x, unpacked_y);
+  if (sum_xx && sum_xy) {
+    *sum_xx += xx;
+    *sum_xy += xy;
+  } else {
+    sum_xx = xx;
+    sum_xy = xy;
+  }
+}
+
+std::vector<std::vector<double>>
+IncrementalLinearLeastSquares::GetBestSolution() const {
+  RTC_CHECK(sum_xx && sum_xy) << "No observations have been added";
+  return ToVectorMatrix(GaussianElimination(ConvertTo<double>(*sum_xx),
+                                            ConvertTo<double>(*sum_xy)));
+}
+
+}  // namespace test
+}  // namespace webrtc