diff --git a/bwe/BandwidthEstimator.cpp b/bwe/BandwidthEstimator.cpp index 944a61433..dbe96040d 100644 --- a/bwe/BandwidthEstimator.cpp +++ b/bwe/BandwidthEstimator.cpp @@ -133,7 +133,10 @@ void BandwidthEstimator::update(uint32_t packetSize, uint64_t transmitTimeNs, ui const auto kalmanGain = crossCovariance * (1.0 / covDelay); _state = predictedMeanState + kalmanGain * (observedDelay - predictedMeanDelay); _covarianceP = statePredictionCovariance - kalmanGain * covDelay * math::transpose(kalmanGain); + math::makeSymmetric(_covarianceP); + assert(math::isSymmetric(_covarianceP)); assert(math::isValid(_covarianceP)); + _observedDelay = observedDelay; _state(0) = std::max(0.0, std::min(_state(0), _config.maxNetworkQueue * 8)); @@ -257,7 +260,8 @@ void BandwidthEstimator::generateSigmaPoints(const math::Matrix& stat const math::Matrix& processNoise, std::array, SIGMA_POINTS>& sigmaPoints) { - auto squareRoot = math::choleskyDecompositionLL(covP); + static const auto seed = covP.I() * 0.0000001; // will make it positive definite + const auto squareRoot = math::choleskyDecompositionLL(covP + seed); sigmaPoints[0] = _state; int startIndex = 1; diff --git a/bwe/Research.md b/bwe/Research.md index aece085c5..bf0f20578 100644 --- a/bwe/Research.md +++ b/bwe/Research.md @@ -1,11 +1,15 @@ -# Tuning the filter +# Tuning the bandwidth estimation filter The Unscented Kalman Filter has a state of current queue length Q, bandwidth BW, and clock offset CO. These three values evolve over time to predict the observed absolute difference between local receive time and remote transmit time of packets. What makes this a bit sensitive is that higher delay can be explained -by lower BW, higher CO or larger Q. If any of process noise values have a bad setting, the state +by lower BW, higher CO or larger Q. If any of process noise values has a bad setting, the state may start to drift in the wrong direction. +The state units are bits, kbps, ms. +The reason for this is to condition the matrix into having comparable magnitudes in the elements. +There are matrix addition, multiplcation and decomposition that will cause + ## Parameter tuning experiments ### Tuning clock drift diff --git a/math/Matrix.h b/math/Matrix.h index 56c1d5609..ec6a40328 100644 --- a/math/Matrix.h +++ b/math/Matrix.h @@ -12,6 +12,8 @@ namespace math template class Matrix { + static_assert(M > 0 && N > 0); + public: Matrix() { @@ -52,7 +54,22 @@ class Matrix return *this; } - bool operator==(const Matrix& other) const { return 0 == std::memcmp(&_m, &other._m, sizeof(_m)); } + bool operator==(const Matrix& other) const + { + for (uint32_t i = 0; i < M; ++i) + { + for (uint32_t j = 0; j < N; ++j) + { + if (_m[i][j] != other._m[i][j]) + { + return false; + } + } + } + return true; + } + + bool operator!=(const Matrix& other) const { return !(*this == other); } Matrix& operator+=(const Matrix& m2) { @@ -78,6 +95,19 @@ class Matrix return *this; } + Matrix operator-() const + { + Matrix n; + for (uint32_t i = 0; i < M; ++i) + { + for (uint32_t j = 0; j < N; ++j) + { + n(i, j) = -_m[i][j]; + } + } + return n; + } + inline T& operator()(uint32_t r, uint32_t c = 0) { assert(r < M && c < N); @@ -99,6 +129,16 @@ class Matrix return v; } + Matrix getRow(uint32_t i) const + { + Matrix v; + for (uint32_t k = 0; k < N; ++k) + { + v(0, k) = _m[i][k]; + } + return v; + } + void setColumn(uint32_t j, const Matrix& v) { for (uint32_t k = 0; k < M; ++k) @@ -110,10 +150,76 @@ class Matrix constexpr int columns() const { return N; } constexpr int rows() const { return M; } + static Matrix I() + { + static_assert(M == N); + Matrix m; + for (auto i = 0u; i < M; ++i) + { + m(i, i) = 1.0; + } + return m; + } + + T absSum() const + { + T s = 0; + for (auto i = 0u; i < M; ++i) + { + for (auto j = 0u; j < M; ++j) + { + s += std::abs(_m[i][j]); + } + } + return s; + } + + void swapRows(uint32_t a, uint32_t b) + { + for (auto i = 0u; i < N; ++i) + { + std::swap(_m[a][i], _m[b][i]); + } + } + + // return which row was swapped, or same rowColumn if not swapped + uint32_t pivot(uint32_t rowColumn) + { + if (rowColumn >= M - 1 || rowColumn >= N - 1) + { + return rowColumn; + } + + T maxElementValue = std::abs(_m[rowColumn][rowColumn]); + uint32_t maxElementRow = rowColumn; + + for (auto k = rowColumn + 1; k < M; k++) + { + auto value = std::abs(_m[k][rowColumn]); + + if (value > maxElementValue) + { + maxElementValue = value; + maxElementRow = k; + } + } + + if (maxElementRow != rowColumn) + { + swapRows(rowColumn, maxElementRow); + return maxElementRow; + } + + return rowColumn; + } + private: T _m[M][N]; }; +template +using Vector = Matrix; + template Matrix operator*(const Matrix& m1, const Matrix& m2) { @@ -222,24 +328,28 @@ Matrix choleskyDecompositionLL(const Matrix& m) { for (uint32_t k = 0; k < j; k++) { - sum += (lowerTriangular(i, k) * lowerTriangular(i, k)); + sum += lowerTriangular(i, k) * lowerTriangular(i, k); } - auto d = m(i, i) - sum; + const auto d = m(i, i) - sum; lowerTriangular(j, j) = std::sqrt(std::abs(d)); } else { for (uint32_t k = 0; k < j; k++) { - sum += (lowerTriangular(i, k) * lowerTriangular(j, k)); + sum += lowerTriangular(i, k) * lowerTriangular(j, k); } - auto d = m(i, j) - sum; - if (lowerTriangular(j, j) == 0) + const auto d = m(i, j) - sum; + if (lowerTriangular(j, j) != 0) { - assert(std::abs(d) < 1E-20); + lowerTriangular(i, j) = d / lowerTriangular(j, j); + } + else + { + assert(std::abs(d) < 1E-10); + lowerTriangular(i, j) = 0; } - lowerTriangular(i, j) = (lowerTriangular(j, j) != 0 ? d / lowerTriangular(j, j) : 0); } } } @@ -282,4 +392,564 @@ bool isValid(const Matrix& m) return true; } +// Frobenius norm of matrix, and "norm" of vector +template +T norm(const Matrix& m) +{ + T n = 0; + for (auto i = 0u; i < M; ++i) + { + for (auto j = 0u; j < N; ++j) + { + n += m(i, j) * m(i, j); + } + } + return std::sqrt(n); +} + +template +T min(const Matrix& m) +{ + T n = m(0, 0); + for (auto i = 0u; i < M; ++i) + { + for (auto j = 0u; j < N; ++j) + { + n = std::min(n, m(i, j)); + } + } + return n; +} + +template +T max(const Matrix& m) +{ + T n = m(0, 0); + for (auto i = 0u; i < M; ++i) + { + for (auto j = 0u; j < N; ++j) + { + n = std::max(n, m(i, j)); + } + } + return n; +} + +template +bool isSymmetric(const Matrix& m) +{ + for (auto i = 1u; i < M; ++i) + { + for (auto j = 0u; j < i; ++j) + { + if (m(i, j) != m(j, i)) + { + return false; + } + } + } + return true; +} + +template +void makeSymmetric(Matrix& m) +{ + for (auto i = 1u; i < M; ++i) + { + for (auto j = 0u; j < i; ++j) + { + if (m(i, j) != m(j, i)) + { + const auto v = (m(i, j) + m(j, i)) / 2.0; + m(i, j) = v; + m(j, i) = v; + } + } + } +} + +// swaps selected row with a lower row to make sure smallest absolute value is at row cr column cr. +// Multiplies the swapped row with -1 to preserve determinant +template +void pivotRow(Matrix& m, uint32_t cr) +{ + static_assert(std::is_signed::value); + static_assert(N >= M); + + auto candidatePos = m.pivot(cr); + if (candidatePos == cr) + { + return; + } + + bool signIsDifferent = std::signbit(m(cr, cr)) != std::signbit(m(candidatePos, cr)); + for (auto j = 0u; j < N; ++j) + { + if (signIsDifferent) + { + m(cr, j) *= -1; + } + else + { + m(candidatePos, j) *= -1; + } + } + + return; +} + +template +void gaussianElimination(Matrix& m) +{ + // along diagonal + for (auto i = 0u; i < M - 1; ++i) + { + pivotRow(m, i); + + if (m(i, i) == 0) + { + continue; + } + + for (auto r = i + 1; r < M; ++r) + { + const auto f = m(r, i) / m(i, i); + m(r, i) = 0; + for (auto j = i + 1; j < N; ++j) + { + m(r, j) -= m(i, j) * f; + } + } + } +} + +template +Matrix eleminateRowColumn(const Matrix& m, uint32_t r, uint32_t c) +{ + Matrix em; + + uint32_t a = 0; + for (auto i = 0u; i < M; ++i) + { + if (i == r) + { + continue; + } + + uint32_t b = 0; + for (auto j = 0u; j < N; ++j) + { + if (j == c) + { + continue; + } + em(a, b) = m(i, j); + ++b; + } + ++a; + } + + return em; +} + +template +T det(const Matrix& m) +{ + return m(0, 0); +} + +template +T det(const Matrix& m) +{ + const auto a = m(0, 0) * m(1, 1); + const auto b = m(0, 1) * m(1, 0); + + return a - b; +} + +// determinant using Gaussian elemination, O(n2) +template +T det(Matrix m) +{ + gaussianElimination(m); + + T a = m(0, 0); + for (auto i = 1u; i < M; ++i) + { + a *= m(i, i); + } + + return a; +} + +template +Matrix principalSubMatrix(const Matrix& m) +{ + static_assert(R < M && R < N); + Matrix psm; + for (auto i = 0u; i < R; ++i) + { + for (auto j = 0u; j < R; ++j) + { + psm(i, j) = m(i, j); + } + } + return psm; +} + +template +T leadingPrincipalMinor(const Matrix& m) +{ + return det(principalSubMatrix(m)); +} + +namespace +{ + +} // namespace + +// m is positive definite if all leading principal minors are positive +// If we do Gaussian elimination first, all principal minors are positive if diagonal elements are positive +template +bool isPositiveDefinite(const Matrix& m) +{ + if (!isSymmetric(m)) + { + return false; + } + + auto pminors = principalMinors(m); + + for (auto i = 0u; i < M; ++i) + { + T v = pminors(i); + if (v <= 0) + { + return false; + } + } + + return true; +} + +// check is positive semi definite by Gaussian elimination +template +bool isPositiveSemiDefinite(const Matrix& m, T tolerance = 1e-14) +{ + if (!isSymmetric(m)) + { + return false; + } + + auto pminors = principalMinors(m); + for (auto i = 0u; i < M; ++i) + { + T v = pminors(i); + if (v < 0 && -v > tolerance) + { + return false; + } + } + + return true; +} + +template +math::Matrix principalMinors(const math::Matrix& m) +{ + return m; +} + +template +math::Matrix principalMinors(const math::Matrix& m) +{ + math::Matrix x; + + auto xs = principalMinors(principalSubMatrix(m)); + + for (auto i = 0u; i < xs.rows(); ++i) + { + x(i) = xs(i); + } + + x(M - 1) = det(m); + + return x; +} + +template +void randomizePositive(Matrix& m, T maxValue) +{ + for (auto i = 0u; i < M; ++i) + { + for (auto j = 0u; j < N; ++j) + { + m(i, j) = maxValue * (rand() % 10000) / 10000.0; + } + } +} + +template +void randomize(Matrix& m, T valueRange) +{ + for (auto i = 0u; i < M; ++i) + { + for (auto j = 0u; j < N; ++j) + { + m(i, j) = valueRange * (rand() % 10000) / 10000.0 - valueRange / 2; + } + } +} + +template +bool isLowerTriangular(const Matrix& m) +{ + for (uint32_t i = 0; i < M; ++i) + { + for (uint32_t j = i + 1; j < M; ++j) + { + if (m(i, j) != 0) + { + return false; + } + } + } + return true; +} + +template +void mulDiagonal(Matrix& m, T scalar) +{ + for (auto i = 0u; i < M; ++i) + { + m(i, i) *= scalar; + } +} + +// Bread and butter for Matrix operations. This output can be used to calculate determinant, invert matrix or solve +// linear equation system. +// Returns (L-E)+U such that P*A=L*U. +template +math::Matrix decomposePLU(math::Matrix m, + math::Matrix& P, + uint32_t& rowSwaps, + const T tolerance = 1e-16) +{ + rowSwaps = 0; + + P = P.I(); + + // iterate rows + for (auto i = 0u; i < M; ++i) + { + auto pivotRow = m.pivot(i); + if (pivotRow != i) + { + P.swapRows(i, pivotRow); + ++rowSwaps; + } + + if (std::abs(m(i, i)) < tolerance) + { + return math::Matrix(); // bad matrix + } + + for (auto j = i + 1; j < M; ++j) + { + if (m(i, i) != 0) + { + m(j, i) /= m(i, i); + } + else + { + m(j, i) = 0; + } + + for (auto k = i + 1; k < M; ++k) + { + m(j, k) -= m(j, i) * m(i, k); + } + } + } + return m; +} + +template +math::Matrix decomposePLU(const math::Matrix& m, math::Matrix& P) +{ + uint32_t dummy; + return decomposePLU(m, P, dummy); +} + +template +void splitPLU(const math::Matrix& lue, math::Matrix& L, math::Matrix& U) +{ + L = math::Matrix(); + U = math::Matrix(); + + for (auto i = 0u; i < M; ++i) + { + L(i, i) = 1.0; + U(i, i) = lue(i, i); + for (auto j = 0u; j < i; ++j) + { + L(i, j) = lue(i, j); + } + for (auto j = i + 1; j < M; ++j) + { + U(i, j) = lue(i, j); + } + } +} + +template +T detLU(const Matrix& m) +{ + Matrix P; + uint32_t swaps = 0; + auto lu = decomposePLU(m, P, swaps); + + T d = lu(0, 0); + for (auto i = 1u; i < M; i++) + d *= lu(i, i); + + return (swaps % 2) == 0 ? d : -d; +} + +template +T detByCoFactors(const Matrix& m) +{ + return det(m); +} + +template +T detByCoFactors(const Matrix& m) +{ + return det(m); +} + +// O(n!) calculation with cofactors and minors along first row +template +T detByCoFactors(const Matrix& m) +{ + T d = 0; + T sgn = 1; // start at (0,0) is (-1)^(1+1) + for (auto i = 0u; i < M; ++i) + { + if (m(0, i) != 0) + { + const auto subDeterminant = detByCoFactors(eleminateRowColumn(m, 0, i)); + const auto p = sgn * m(0, i) * subDeterminant; + + if (std::signbit(p) != std::signbit(d) && (std::nextafter(d, -p) == p || std::nextafter(-p, d) == d)) + { + d = 0; + } + else + { + d += p; + } + } + sgn = sgn > 0 ? -1. : 1.; + } + return d; +} + +template +math::Matrix invertLU(const math::Matrix& LUe, const math::Matrix& P) +{ + math::Matrix iMatrix = P; + + for (auto col = 0u; col < N; ++col) + { + for (auto i = 0u; i < N; ++i) + { + for (auto k = 0u; k < i; ++k) + { + iMatrix(i, col) -= LUe(i, k) * iMatrix(k, col); + } + } + + for (auto i = N - 1; i != (0u - 1); --i) + { + for (auto k = i + 1; k < N; ++k) + { + iMatrix(i, col) -= LUe(i, k) * iMatrix(k, col); + } + + iMatrix(i, col) /= LUe(i, i); + } + } + return iMatrix; +} + +template +math::Matrix solveLU(const math::Matrix& lue, + const math::Matrix& P, + const math::Matrix& b) +{ + math::Matrix x; + auto pb = P * b; + for (auto i = 0u; i < M; ++i) + { + x(i) = pb(i); + + for (auto k = 0u; k < i; ++k) + { + x(i) -= lue(i, k) * x(k); + } + } + + for (auto i = 1u; i <= M; ++i) + { + for (auto k = M - i + 1; k < M; ++k) + { + x(M - i) -= lue(M - i, k) * x(k); + } + x(M - i) /= lue(M - i, M - i); + } + + return x; +} + +template +math::Matrix solve(math::Matrix a, const math::Matrix& b) +{ + auto ab = augment(a, b); + gaussianElimination(ab); + + math::Matrix x; + for (auto k = 1u; k <= M; ++k) + { + const auto i = M - k; + T v = ab(i, M); + for (auto j = i + 1; j < M; ++j) + { + v -= ab(i, j) * x(j); + } + x(i) = v / ab(i, i); + } + + return x; +} + +template +math::Matrix augment(const math::Matrix& m, const math::Matrix& b) +{ + math::Matrix a; + for (auto i = 0u; i < M; ++i) + { + for (auto j = 0u; j < N; ++j) + { + a(i, j) = m(i, j); + } + + for (auto j = 0u; j < W; ++j) + { + a(i, j + N) = b(i, j); + } + } + + return a; +} } // namespace math diff --git a/rtp/JitterEstimator.h b/rtp/JitterEstimator.h index 31b106a9d..ae46fa8f1 100644 --- a/rtp/JitterEstimator.h +++ b/rtp/JitterEstimator.h @@ -1,5 +1,4 @@ #pragma once -#include "math/Matrix.h" #include "math/WelfordVariance.h" #include "rtp/RtpDelayTracker.h" #include "utils/Trackers.h" diff --git a/test/bwe/EstimatorReRun.cpp b/test/bwe/EstimatorReRun.cpp index 914f4bf48..67392e502 100644 --- a/test/bwe/EstimatorReRun.cpp +++ b/test/bwe/EstimatorReRun.cpp @@ -144,16 +144,16 @@ uint32_t identifyAudioSsrc(logger::PacketLogReader& reader) } } // namespace -class BweRerunTrace : public testing::TestWithParam +class BweRerun : public testing::TestWithParam { }; -TEST_P(BweRerunTrace, fromTrace) +TEST_P(BweRerun, DISABLED_fromTrace) { bwe::Config config; config.congestion.cap.ratio = 0.5; - std::string trace = GetParam(); + auto trace = GetParam(); if (trace.empty()) { @@ -237,33 +237,33 @@ TEST_P(BweRerunTrace, fromTrace) logger::info("%" PRIu64 "s finished at %.3fkbps", "", (item.receiveTimestamp - start) / utils::Time::sec, prevBw); } -INSTANTIATE_TEST_SUITE_P(DISABLED_BweRerun, - BweRerunTrace, - testing::Values("Transport-39_tcp", - "Transport-58_tcp", - "Transport-86_tcp_1ploss", +INSTANTIATE_TEST_SUITE_P(BweReruns, + BweRerun, + testing::Values("Transport-4-wifi", "Transport-105_tcp_1ploss", - "Transport-22-4G-2.3Mbps", - "Transport-1094-4G", + "Transport-22-3G", + "Transport-30_oka", "Transport-3644-wifi", + "Transport-39_tcp", + "Transport-4735-4G", + "Transport-48_80_3G", "Transport-62-4G", - "Transport-3887-wifi", - "Transport-3629", - "Transport-14-wifi", - "Transport-4-wifi", - "Transport-6-4G-1-5Mbps", - "Transport-30-3G-1Mbps", + "Transport-1094-4G", + "Transport-22-4G-2.3Mbps", "Transport-32_Idre", - "Transport-30_oka", + "Transport-37", "Transport-48_50_3G", + "Transport-58_tcp", + "Transport-86_tcp_1ploss", + "Transport-14-wifi", + "Transport-30-3G-1Mbps", + "Transport-3629", + "Transport-3887-wifi", + "Transport-44-clkdrift", "Transport-48_60_3G", - "Transport-48_80_3G", - "Transport-22-3G", - "Transport-62-4G", - "Transport-4735-4G", - "Transport-44-clkdrift")); + "Transport-6-4G-1-5Mbps")); -class BweRerunLimit : public testing::TestWithParam +class BweRerunLimit : public testing::TestWithParam> { }; @@ -271,116 +271,137 @@ TEST_P(BweRerunLimit, DISABLED_limitedLink) { bwe::Config config; - std::array trace = {"Transport-5", "Transport-20", "Transport-17"}; + auto trace = std::get<0>(GetParam()); + auto linkCapacity = std::get<1>(GetParam()); + memory::PacketPoolAllocator allocator(8092, "rerun"); - fakenet::NetworkLink link("EstimatorReRunLink", GetParam(), 1950 * 1024, 3000); + fakenet::NetworkLink link("EstimatorReRunLink", linkCapacity, 1950 * 1024, 3000); link.setLossRate(0); uint64_t wallClock = 0; const char* formatLine = "%u, %.1f,%.1f,%.f,%.3f,%.3f,%u,%u,%.6f, %.2f"; const char* legend = "time, bwo, bw, q, co, delay, size, seqno, stime, rate"; - for (size_t t = 0; t < trace.size() && !trace[t].empty(); ++t) - { - bwe::BandwidthEstimator estimator(config); - rtp::SendTimeDial sendTimeDial; - std::string path = "./_bwelogs"; - logger::PacketLogReader reader(::fopen((path + "/" + trace[t]).c_str(), "r")); - logger::PacketLogItem item; + bwe::BandwidthEstimator estimator(config); + rtp::SendTimeDial sendTimeDial; + std::string path = "./_bwelogs"; + logger::PacketLogReader reader(::fopen((path + "/" + trace).c_str(), "r")); + logger::PacketLogItem item; - CsvWriter csvLog(("./_ssdatall/" + trace[t] + std::to_string(GetParam()) + ".csv").c_str()); + CsvWriter csvLog(("./_ssdata/" + trace + std::to_string(linkCapacity) + "All.csv").c_str()); - csvLog.writeLine("%s", legend); + csvLog.writeLine("%s", legend); - logger::info("processing %s", "bweRerun", trace[t].c_str()); - for (int i = 0; reader.getNext(item); ++i) + logger::info("processing %s", "bweRerun", trace.c_str()); + for (int i = 0; reader.getNext(item); ++i) + { + if (i == 0) { - if (i == 0) - { - wallClock = sendTimeDial.toAbsoluteTime(item.transmitTimestamp, utils::Time::minute); - auto packet = memory::makeUniquePacket(allocator, &item, sizeof(item)); - packet->setLength(item.size); + wallClock = sendTimeDial.toAbsoluteTime(item.transmitTimestamp, utils::Time::minute); + auto packet = memory::makeUniquePacket(allocator, &item, sizeof(item)); + packet->setLength(item.size); - link.push(std::move(packet), wallClock); - continue; - } + link.push(std::move(packet), wallClock); + continue; + } - auto sendTime = sendTimeDial.toAbsoluteTime(item.transmitTimestamp, wallClock + utils::Time::sec * 10); - assert(sendTime >= wallClock); - for (;;) + auto sendTime = sendTimeDial.toAbsoluteTime(item.transmitTimestamp, wallClock + utils::Time::sec * 10); + assert(sendTime >= wallClock); + for (;;) + { + auto minAdvance = std::min(utils::Time::diff(wallClock, sendTime), link.timeToRelease(wallClock)); + assert(minAdvance >= 0); + wallClock += minAdvance; + for (auto packet = link.pop(wallClock); packet; packet = link.pop(wallClock)) { - auto minAdvance = std::min(utils::Time::diff(wallClock, sendTime), link.timeToRelease(wallClock)); - assert(minAdvance >= 0); - wallClock += minAdvance; - for (auto packet = link.pop(wallClock); packet; packet = link.pop(wallClock)) - { - if (packet->getLength() < 1500) - { - auto* packetItem = reinterpret_cast(packet->get()); - - auto packetSendTime = sendTimeDial.toAbsoluteTime(packetItem->transmitTimestamp, - wallClock + utils::Time::sec * 10); - estimator.update(packetItem->size, packetSendTime, wallClock); - - auto bw = estimator.getEstimate(item.receiveTimestamp); - auto state = estimator.getState(); - auto delay = estimator.getDelay(); - csvLog.writeLine(formatLine, - wallClock / utils::Time::ms, - bw, - state(1), - state(0) / 8, - state(2), - delay, - packetItem->size, - packetItem->sequenceNumber, - double(item.transmitTimestamp >> 18) + double(item.transmitTimestamp & 0x3FFFF) / 262144, - std::min(9000.0, estimator.getReceiveRate(wallClock))); - } - } - if (wallClock == sendTime) + if (packet->getLength() < 1500) { - auto packet = memory::makeUniquePacket(allocator, &item, sizeof(item)); - packet->setLength(item.size); - link.push(std::move(packet), wallClock); - - break; + auto* packetItem = reinterpret_cast(packet->get()); + + auto packetSendTime = + sendTimeDial.toAbsoluteTime(packetItem->transmitTimestamp, wallClock + utils::Time::sec * 10); + estimator.update(packetItem->size, packetSendTime, wallClock); + + auto bw = estimator.getEstimate(item.receiveTimestamp); + auto state = estimator.getState(); + auto delay = estimator.getDelay(); + csvLog.writeLine(formatLine, + wallClock / utils::Time::ms, + bw, + state(1), + state(0) / 8, + state(2), + delay, + packetItem->size, + packetItem->sequenceNumber, + double(item.transmitTimestamp >> 18) + double(item.transmitTimestamp & 0x3FFFF) / 262144, + std::min(9000.0, estimator.getReceiveRate(wallClock))); } } - } - - double lastEstimate = estimator.getEstimate(wallClock); - for (wallClock += link.timeToRelease(wallClock);; wallClock += link.timeToRelease(wallClock)) - { - auto packet = link.pop(wallClock); - if (!packet) + if (wallClock == sendTime) { + auto packet = memory::makeUniquePacket(allocator, &item, sizeof(item)); + packet->setLength(item.size); + link.push(std::move(packet), wallClock); + break; } - auto* packetItem = reinterpret_cast(packet->get()); - auto sendTime = - sendTimeDial.toAbsoluteTime(packetItem->transmitTimestamp, wallClock + utils::Time::sec * 10); - estimator.update(packetItem->size, sendTime, wallClock); - auto bw = estimator.getEstimate(item.receiveTimestamp); - auto state = estimator.getState(); - auto delay = estimator.getDelay(); - csvLog.writeLine(formatLine, - wallClock / utils::Time::ms, - bw, - state(1), - state(0) / 8, - state(2), - delay, - packetItem->size, - packetItem->sequenceNumber, - double(item.transmitTimestamp >> 18) + double(item.transmitTimestamp & 0x3FFFF) / 262144, - std::min(9000.0, estimator.getReceiveRate(wallClock))); + } + } - lastEstimate = estimator.getEstimate(wallClock); + double lastEstimate = estimator.getEstimate(wallClock); + for (wallClock += link.timeToRelease(wallClock);; wallClock += link.timeToRelease(wallClock)) + { + auto packet = link.pop(wallClock); + if (!packet) + { + break; } - logger::info("finished at %.3fkbps", "", lastEstimate); + auto* packetItem = reinterpret_cast(packet->get()); + auto sendTime = sendTimeDial.toAbsoluteTime(packetItem->transmitTimestamp, wallClock + utils::Time::sec * 10); + estimator.update(packetItem->size, sendTime, wallClock); + auto bw = estimator.getEstimate(item.receiveTimestamp); + auto state = estimator.getState(); + auto delay = estimator.getDelay(); + csvLog.writeLine(formatLine, + wallClock / utils::Time::ms, + bw, + state(1), + state(0) / 8, + state(2), + delay, + packetItem->size, + packetItem->sequenceNumber, + double(item.transmitTimestamp >> 18) + double(item.transmitTimestamp & 0x3FFFF) / 262144, + std::min(9000.0, estimator.getReceiveRate(wallClock))); + + lastEstimate = estimator.getEstimate(wallClock); } + logger::info("finished at %.3fkbps", "", lastEstimate); } INSTANTIATE_TEST_SUITE_P(BwessRerun, BweRerunLimit, - testing::Values(4500, 4800, 4900, 5000, 5100, 5200, 5600, 6500, 7000, 10000, 20000)); + testing::Combine(testing::Values("Transport-4-wifi", + "Transport-105_tcp_1ploss", + "Transport-22-3G", + "Transport-30_oka", + "Transport-3644-wifi", + "Transport-39_tcp", + "Transport-4735-4G", + "Transport-48_80_3G", + "Transport-62-4G", + "Transport-1094-4G", + "Transport-22-4G-2.3Mbps", + "Transport-32_Idre", + "Transport-37", + "Transport-48_50_3G", + "Transport-58_tcp", + "Transport-86_tcp_1ploss", + "Transport-14-wifi", + "Transport-30-3G-1Mbps", + "Transport-3629", + "Transport-3887-wifi", + "Transport-44-clkdrift", + "Transport-48_60_3G", + "Transport-6-4G-1-5Mbps"), + testing::Values(4500, 4800, 4900, 5000, 5100, 5200, 5600, 6500, 7000, 10000, 20000))); diff --git a/test/bwe/MatrixTests.cpp b/test/bwe/MatrixTests.cpp index 2c8951c0e..5d821ac2b 100644 --- a/test/bwe/MatrixTests.cpp +++ b/test/bwe/MatrixTests.cpp @@ -1,6 +1,8 @@ +#include "logger/Logger.h" #include "math/Matrix.h" #include #include + using namespace math; TEST(MatrixTest, mul) { @@ -18,7 +20,7 @@ TEST(MatrixTest, mul) auto m3 = m * transpose(m2); EXPECT_EQ(m3.columns(), m2.rows()); - Matrix v({5.2, 3.5, 2.6, 1.2}); + Matrix v({5.2, 3.5, 2.6, 1.2}); auto m4 = v * transpose(v); EXPECT_EQ(m4.columns(), m4.rows()); EXPECT_EQ(m4.columns(), v.rows()); @@ -26,12 +28,20 @@ TEST(MatrixTest, mul) auto mk = choleskyDecompositionLL(m4); EXPECT_EQ(mk(0, 0), sqrt(m4(0, 0))); } - + TEST(MatrixTest, cholesky) { math::Matrix v({5.6, -156000.0, 19.9}); - auto m = v * math::transpose(v); + auto m = math::outerProduct(v); + + EXPECT_NEAR(det(math::principalSubMatrix(m)), 0.0, 0.0002); + auto a0 = det(m); + auto c0 = detLU(m); + EXPECT_NEAR(c0, 0.0, 6e-17); + EXPECT_NEAR(a0, 0.0, 6e-17); + // should be psd, but representaion errors causes det to be slightly negative + EXPECT_TRUE(math::isPositiveSemiDefinite(m, 0.001)); auto L = choleskyDecompositionLL(m); EXPECT_EQ(L.getColumn(0), v); auto c = L * transpose(L); @@ -46,9 +56,302 @@ TEST(MatrixTest, choleskyStability) {5.5034655381268649E-19, 102615.51002615986, -3.4396659613292905E-20}, {-3.9095946203353344E-29, -3.4396659613292905E-20, 2.443496637709584E-30}}); + ASSERT_TRUE(math::isSymmetric(m)); + EXPECT_TRUE(isPositiveSemiDefinite(m)); + + ASSERT_GE(math::det(m), 0); auto L = choleskyDecompositionLL(m); auto c = L * transpose(L); EXPECT_FALSE(m == c); + EXPECT_LT(math::norm(m - c), 1.E-10); EXPECT_TRUE(math::isValid(L)); // should check determinant of m - c -} \ No newline at end of file +} + +TEST(MatrixTest, determinant) +{ + math::Matrix m; + + const auto acceptableRelativeError = 0.00000001; + for (auto k = 0u; k < 15000; ++k) + { + randomize(m, 1000.0); + + const auto a0 = detByCoFactors(m); + const auto a1 = detByCoFactors(transpose(m)); + const auto b0 = det(m); + const auto b1 = det(transpose(m)); + const auto d1 = detLU(m); + const auto d2 = detLU(transpose(m)); + + const auto absError = acceptableRelativeError * std::abs(d1 + d2) / 2; + EXPECT_NEAR(d1, d2, absError); + EXPECT_NEAR(b0, b1, absError); + EXPECT_NEAR(a0, a1, absError); + } +} + +TEST(MatrixTest, submatrix) +{ + math::Matrix m; + math::randomize(m, 1000.0); + + auto s2 = math::principalSubMatrix(m); + EXPECT_EQ(m(0, 0), s2(0, 0)); + EXPECT_EQ(m(0, 1), s2(0, 1)); + EXPECT_EQ(m(1, 0), s2(1, 0)); + EXPECT_EQ(m(1, 1), s2(1, 1)); +} + +TEST(MatrixTest, posDefininte) +{ + math::Matrix m2({{10.0, 5., 2.}, {5., 3., 2.}, {2., 2., 3.}}); + ASSERT_TRUE(isSymmetric(m2)); + EXPECT_TRUE(isPositiveDefinite(m2)); +} + +TEST(MatrixTest, cholesky3) +{ + math::Matrix v({10, 5, 3}); + auto m = math::outerProduct(v); + + ASSERT_TRUE(math::isSymmetric(m)); + EXPECT_TRUE(math::isPositiveSemiDefinite(m)); + + auto L = choleskyDecompositionLL(m); + auto c = L * transpose(L); + EXPECT_EQ(m, c); +} + +TEST(MatrixTest, choleskyOuterProduct) +{ + const auto seed = math::Matrix::I() * 0.0000001; + + for (auto i = 0u; i < 10000; ++i) + { + math::Matrix v({10, 5, 3}); + math::randomize(v, 10000.0); + auto m = math::outerProduct(v); + + EXPECT_TRUE(isPositiveDefinite(m + seed)); + + const auto L = choleskyDecompositionLL(m); + EXPECT_TRUE(math::isLowerTriangular(L)); + EXPECT_EQ(m, L * transpose(L)); + } +} + +TEST(MatrixTest, PDtest) +{ + const auto seed = math::Matrix::I() * 0.0000001; + double data[][3] = {// + {-2418.000000, 2019.000000, 471.000000}, + {2778.000000, 4342.000000, 4969.000000}, + {-4223.000000, 2884.000000, 2390.000000}, + {3809.000000, -2757.000000, -1085.000000}, + {2697.000000, -1958.000000, -3622.000000}, + {2514.000000, 3210.000000, -3374.000000}, + {3206.000000, 4751.000000, 2886.000000}, + {-1684.000000, -2378.000000, -4991.000000}, + {3262.000000, 4944.000000, -1311.000000}, + {2779.000000, -2046.000000, 1061.000000}, + {-3189.000000, -4757.000000, -1653.000000}, + {3566.000000, -2718.000000, -2750.000000}, + {-3384.000000, -4629.000000, -4016.000000}, + {3425.000000, 4602.000000, 2874.000000}, + {-1459.000000, -2307.000000, -3262.000000}, + {-3847.000000, 2437.000000, -3798.000000}, + {-2956.000000, 1891.000000, 522.000000}, + {2725.000000, -2023.000000, -3784.000000}, + {3499.000000, -2884.000000, -1906.000000}, + {-2841.000000, -4542.000000, 593.000000}, + {-2549.000000, 1967.000000, 752.000000}, + {-3158.000000, -4790.000000, -4507.000000}, + {1907.000000, -1366.000000, 2179.000000}, + {-903.000000, -1092.000000, 4329.000000}, + {-1388.000000, 983.000000, -1764.000000}, + {2166.000000, 3134.000000, 1795.000000}, + {-3913.000000, 2419.000000, -3397.000000}, + {2609.000000, 4471.000000, 2350.000000}, + {-2440.000000, -3375.000000, 4311.000000}, + {2928.000000, 4678.000000, -4407.000000}, + {1193.000000, 1756.000000, -1151.000000}, + {-3619.000000, -4628.000000, 3646.000000}, + {-3212.000000, -4704.000000, 4656.000000}, + {-1650.000000, -2381.000000, 4462.000000}, + {-2269.000000, -3327.000000, 2056.000000}, + {3430.000000, 4589.000000, -4473.000000}, + {-1006.000000, 660.000000, 1035.000000}, + {-3533.000000, -4649.000000, 1754.000000}, + {3747.000000, -2536.000000, 3672.000000}, + {-2988.000000, -4716.000000, -2542.000000}, + {3294.000000, 4731.000000, -4378.000000}}; + + for (auto rawVector : data) + { + math::Matrix v({rawVector[0], rawVector[1], rawVector[2]}); + auto m = math::outerProduct(v); + + EXPECT_TRUE(isPositiveSemiDefinite(m)); + EXPECT_TRUE(isPositiveDefinite(m + seed)); + + const auto L = choleskyDecompositionLL(m); + EXPECT_TRUE(math::isLowerTriangular(L)); + EXPECT_EQ(m, L * transpose(L)); + } +} + +TEST(MatrixTest, principalDeterminants) +{ + math::Matrix v({337, 2614, -2297}); + auto m = math::outerProduct(v); + + auto dx = principalMinors(m); + EXPECT_EQ(dx(0), 337 * 337.0); + EXPECT_EQ(dx(1), 0.0); + EXPECT_EQ(dx(2), 0.0); + EXPECT_TRUE(math::isPositiveSemiDefinite(m)); +} + +TEST(MatrixTest, choleskyOnSum) +{ + for (int i = 0; i < 500; ++i) + { + math::Matrix v({10, 5, 3}); + math::randomize(v, 10000.0); + auto m = math::outerProduct(v); + + EXPECT_TRUE(math::isPositiveSemiDefinite(m)); + EXPECT_FALSE(math::isPositiveDefinite(m)); + + auto sum = m; + for (int j = 0; j < 15; ++j) + { + math::Matrix w; + math::randomize(w, 10.0); + auto smallDiff = math::outerProduct(w); + sum += smallDiff; + } + mulDiagonal(sum, 1.000001); + EXPECT_TRUE(math::isPositiveDefinite(sum)); + EXPECT_TRUE(math::isPositiveSemiDefinite(sum)); + + auto L = math::choleskyDecompositionLL(sum); + EXPECT_TRUE(math::isLowerTriangular(L)); + EXPECT_TRUE(math::isSymmetric(sum)); + if (sum != L * transpose(L)) + { + const auto d = sum - L * transpose(L); + const auto nrm = math::norm(d); + EXPECT_LT(nrm, 1.E-8); + } + } +} + +// Assume two symetric matrices A' and B' created from outer products of vector A and B +// Both A' and B are positive semi definite +// But (A' + B') is positive definite as long as vector elements a1*b2 != a2*b1 +// det(A'+ B') = (a1*b2 - b1*a2)^2, which is always positive as long as inequality above is true +TEST(MatrixTest, semiDefSum2x2) +{ + { + math::Matrix a({4, 3}); + math::Matrix b({8, 6}); + + auto h = math::outerProduct(a) + math::outerProduct(b); + EXPECT_FALSE(math::isPositiveDefinite(h)); + } + + { + math::Matrix a({4, 3}); + math::Matrix b({7, 5}); + + auto h = math::outerProduct(a) + math::outerProduct(b); + EXPECT_TRUE(math::isPositiveDefinite(h)); + } +} + +// For sum of two symmetric 3x3 matrices, the sum is positive semi definite always +// if matrices are produced by outer product +TEST(MatrixTest, semiDefSum3x3) +{ + math::Matrix a({4, 3, 1}); + math::Matrix b({8, 6, 8}); + + auto h = math::outerProduct(a) + math::outerProduct(b); + EXPECT_TRUE(math::isPositiveSemiDefinite(h)); +} + +TEST(MatrixTest, semiDefSum4x4) +{ + math::Matrix a({4, 3, 1, 13}); + math::Matrix b({81, 6, 8, 6}); + + auto h = math::outerProduct(a) + math::outerProduct(b); + auto d = det(h); + EXPECT_TRUE(math::isSymmetric(h)); + EXPECT_NEAR(d, 0, 1e-20); + auto L = math::choleskyDecompositionLL(h); + auto ll = L * transpose(L); + EXPECT_LT(math::norm(h - ll), 1e-7); +} + +TEST(MatrixTest, semiDefSumHuge) +{ + math::Matrix a({4, 3, 1, 13, 56, 1, 2, 5, 4, 2, 4, 8, 3, 6}); + math::Matrix b({81, 6, 8, 6, 3, 8, 1, 9, 4, 2, 6, 23, 6, 4}); + + auto h = math::outerProduct(a) + math::outerProduct(b); + auto d = det(h); + EXPECT_TRUE(math::isSymmetric(h)); + EXPECT_EQ(d, 0); + EXPECT_TRUE(math::isPositiveDefinite(h + h.I() * 0.00001)); + auto L = math::choleskyDecompositionLL(h); + auto ll = L * transpose(L); + EXPECT_LT(math::norm(h - ll), 1e-7); +} + +TEST(MatrixTest, inversion) +{ + math::Matrix a({4, 3, 1, 13}); + auto A = math::outerProduct(a); + A += A.I() * 0.000001; + + auto P = A.I(); + auto lue = math::decomposePLU(A, P); + + math::Matrix L; + math::Matrix U; + math::Matrix E; + + math::splitPLU(lue, L, U); + + EXPECT_LT(math::norm(P * A - L * U), 0.000001); + + auto iA = math::invertLU(lue, P); + + auto identity = A * iA; + auto idiff = identity - identity.I(); + EXPECT_LT(math::norm(idiff), 0.00000001); +} + +TEST(MatrixTest, solveEq) +{ + math::Matrix A({{1, 1, 1}, {0, 2, 5}, {2, 5, -1}}); + math::Matrix b({6, -4, 27}); + + auto ab = augment(A, b); + math::gaussianElimination(ab); + auto x1 = math::solve(A, b); + math::Matrix P; + auto lu = math::decomposePLU(A, P); + auto x = math::solveLU(lu, P, b); + math::Matrix expectedX({5, 3, -2}); + EXPECT_EQ(x, expectedX); + EXPECT_EQ(x1, expectedX); + + auto ia = math::invertLU(lu, P); + auto x2 = ia * b; + EXPECT_LT(math::norm(x2 - expectedX), 0.000000000001); + // auto x = math::solve(lue, P, b); +}