Use 128 bit math for InterpolateOffset
We were being lazy with InterpolateOffset and using the first point as
the integer part of the result, and putting all of the error in the
double. That gets imprecise when we extrapolate long distances, which
is common at the start of log files.
Instead, use the exact integer math, and use the double to capture the
remainder. This is significantly more precise, since we should never
see values over 1.0 in the double.
Change-Id: I02bf511239cd20dbff7d17327db2f665c6aa41d3
Signed-off-by: Austin Schuh <austin.linux@gmail.com>
diff --git a/aos/events/logging/boot_timestamp.h b/aos/events/logging/boot_timestamp.h
index d0fde73..e818fb9 100644
--- a/aos/events/logging/boot_timestamp.h
+++ b/aos/events/logging/boot_timestamp.h
@@ -4,6 +4,7 @@
#include <iostream>
#include "aos/time/time.h"
+#include "glog/logging.h"
namespace aos::logger {
@@ -24,6 +25,18 @@
return {boot, duration - d};
}
+ BootDuration operator-(BootDuration d) const {
+ CHECK_EQ(d.boot, boot);
+ return {boot, duration - d.duration};
+ }
+
+ BootDuration operator+(BootDuration d) const {
+ CHECK_EQ(d.boot, boot);
+ return {boot, duration + d.duration};
+ }
+
+ BootDuration operator/(int x) const { return {boot, duration / x}; }
+
bool operator==(const BootDuration &m2) const {
return boot == m2.boot && duration == m2.duration;
}
diff --git a/aos/network/timestamp_filter.cc b/aos/network/timestamp_filter.cc
index bdfe67f..4ec9e0a 100644
--- a/aos/network/timestamp_filter.cc
+++ b/aos/network/timestamp_filter.cc
@@ -801,24 +801,12 @@
chrono::nanoseconds NoncausalTimestampFilter::InterpolateOffset(
std::tuple<monotonic_clock::time_point, chrono::nanoseconds> p0,
- std::tuple<monotonic_clock::time_point, chrono::nanoseconds> /*p1*/,
- monotonic_clock::time_point /*ta_base*/, double /*ta*/) {
- // For the double variant, we want to split the result up into a large integer
- // portion, and the rest. We want to do this without introducing numerical
- // precision problems.
- //
- // One way would be to carefully compute the integer portion, and then compute
- // the double portion in such a way that the two are guaranteed to add up
- // correctly.
- //
- // The simpler way is to simply just use the offset from p0 as the integer
- // portion, and make the rest be the double portion. It will get us most of
- // the way there for a lot less work, and we can revisit if this breaks down.
- //
- // oa = p0.o + (ta - p0.t) * (p1.o - p0.o) / (p1.t - p0.t)
- // ^^^^
- // TODO(austin): Use 128 bit math and the remainder to be more accurate here.
- return std::get<1>(p0);
+ std::tuple<monotonic_clock::time_point, chrono::nanoseconds> p1,
+ monotonic_clock::time_point ta_base, double ta) {
+ DCHECK_GE(ta, 0.0);
+ DCHECK_LT(ta, 1.0);
+
+ return InterpolateOffset(p0, p1, ta_base);
}
double NoncausalTimestampFilter::InterpolateOffsetRemainder(
@@ -827,15 +815,19 @@
monotonic_clock::time_point ta_base, double ta) {
const chrono::nanoseconds time_in = ta_base - std::get<0>(p0);
const chrono::nanoseconds dt = std::get<0>(p1) - std::get<0>(p0);
+ const chrono::nanoseconds doffset = std::get<1>(p1) - std::get<1>(p0);
- // The remainder then is the rest of the equation.
- //
- // oa = p0.o + (ta - p0.t) * (p1.o - p0.o) / (p1.t - p0.t)
- // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
- // TODO(austin): Use 128 bit math and the remainder to be more accurate here.
- return static_cast<double>(ta + time_in.count()) /
- static_cast<double>(dt.count()) *
- (std::get<1>(p1) - std::get<1>(p0)).count();
+ // Compute the remainder of the division in InterpolateOffset above, and then
+ // use double math to compute it accurately.
+ absl::int128 numerator =
+ absl::int128(time_in.count()) * absl::int128(doffset.count());
+ numerator += numerator > 0 ? absl::int128(dt.count() / 2)
+ : -absl::int128(dt.count() / 2);
+ return static_cast<double>(numerator % absl::int128(dt.count())) /
+ dt.count() +
+ (numerator > 0 ? -0.5 : 0.5) +
+ ta * static_cast<double>(doffset.count()) /
+ static_cast<double>(dt.count());
}
chrono::nanoseconds NoncausalTimestampFilter::BoundOffset(
diff --git a/aos/network/timestamp_filter_test.cc b/aos/network/timestamp_filter_test.cc
index 67afc6e..1669e62 100644
--- a/aos/network/timestamp_filter_test.cc
+++ b/aos/network/timestamp_filter_test.cc
@@ -797,11 +797,9 @@
const monotonic_clock::time_point t1 = e + chrono::nanoseconds(10000);
const chrono::nanoseconds o1 = chrono::nanoseconds(100);
- const double o1d = static_cast<double>(o1.count());
const monotonic_clock::time_point t2 = t1 + chrono::nanoseconds(1000);
const chrono::nanoseconds o2 = chrono::nanoseconds(150);
- const double o2d = static_cast<double>(o2.count());
EXPECT_EQ(NoncausalTimestampFilter::InterpolateOffset(
std::make_tuple(t1, o1), std::make_tuple(t2, o2), t1),
@@ -818,10 +816,10 @@
o2);
EXPECT_EQ(NoncausalTimestampFilter::InterpolateOffset(
std::make_tuple(t1, o1), std::make_tuple(t2, o2), t2, 0.0),
- o1);
+ o2);
EXPECT_EQ(NoncausalTimestampFilter::InterpolateOffsetRemainder(
std::make_tuple(t1, o1), std::make_tuple(t2, o2), t2, 0.0),
- o2d - o1d);
+ 0.0);
EXPECT_EQ(NoncausalTimestampFilter::InterpolateOffset(
std::make_tuple(t1, o1), std::make_tuple(t2, o2),
@@ -830,55 +828,101 @@
EXPECT_EQ(NoncausalTimestampFilter::InterpolateOffset(
std::make_tuple(t1, o1), std::make_tuple(t2, o2),
t1 + chrono::nanoseconds(500), 0.0),
- o1);
+ o1 + chrono::nanoseconds(25));
EXPECT_EQ(NoncausalTimestampFilter::InterpolateOffsetRemainder(
std::make_tuple(t1, o1), std::make_tuple(t2, o2),
t1 + chrono::nanoseconds(500), 0.0),
- 25.);
+ 0.0);
EXPECT_EQ(NoncausalTimestampFilter::InterpolateOffset(
std::make_tuple(t1, o1), std::make_tuple(t2, o2),
t1 + chrono::nanoseconds(-200)),
chrono::nanoseconds(90));
EXPECT_EQ(NoncausalTimestampFilter::InterpolateOffset(
- std::make_tuple(t1, o1), std::make_tuple(t2, o2), t1, -200.),
- o1);
+ std::make_tuple(t1, o1), std::make_tuple(t2, o2),
+ t1 - chrono::nanoseconds(200), 0.0),
+ o1 - chrono::nanoseconds(10));
EXPECT_EQ(NoncausalTimestampFilter::InterpolateOffsetRemainder(
- std::make_tuple(t1, o1), std::make_tuple(t2, o2), t1, -200.),
- -10.);
+ std::make_tuple(t1, o1), std::make_tuple(t2, o2),
+ t1 - chrono::nanoseconds(200), 0.0),
+ 0.0);
EXPECT_EQ(NoncausalTimestampFilter::InterpolateOffset(
std::make_tuple(t1, o1), std::make_tuple(t2, o2),
t1 + chrono::nanoseconds(200)),
chrono::nanoseconds(110));
EXPECT_EQ(NoncausalTimestampFilter::InterpolateOffset(
- std::make_tuple(t1, o1), std::make_tuple(t2, o2), t1, 200.),
- o1);
+ std::make_tuple(t1, o1), std::make_tuple(t2, o2),
+ t1 + chrono::nanoseconds(200), 0.0),
+ o1 + chrono::nanoseconds(10));
EXPECT_EQ(NoncausalTimestampFilter::InterpolateOffsetRemainder(
- std::make_tuple(t1, o1), std::make_tuple(t2, o2), t1, 200.),
- 10.);
+ std::make_tuple(t1, o1), std::make_tuple(t2, o2),
+ t1 + chrono::nanoseconds(200), 0.0),
+ 0.0);
EXPECT_EQ(NoncausalTimestampFilter::InterpolateOffset(
std::make_tuple(t1, o1), std::make_tuple(t2, o2),
t1 + chrono::nanoseconds(800)),
chrono::nanoseconds(140));
EXPECT_EQ(NoncausalTimestampFilter::InterpolateOffset(
- std::make_tuple(t1, o1), std::make_tuple(t2, o2), t1, 800.),
- o1);
+ std::make_tuple(t1, o1), std::make_tuple(t2, o2),
+ t1 + chrono::nanoseconds(800), 0.0),
+ o1 + chrono::nanoseconds(40));
EXPECT_EQ(NoncausalTimestampFilter::InterpolateOffsetRemainder(
- std::make_tuple(t1, o1), std::make_tuple(t2, o2), t1, 800.),
- 40.);
+ std::make_tuple(t1, o1), std::make_tuple(t2, o2),
+ t1 + chrono::nanoseconds(800), 0.0),
+ 0.0);
EXPECT_EQ(NoncausalTimestampFilter::InterpolateOffset(
std::make_tuple(t1, o1), std::make_tuple(t2, o2),
t1 + chrono::nanoseconds(1200)),
chrono::nanoseconds(160));
EXPECT_EQ(NoncausalTimestampFilter::InterpolateOffset(
- std::make_tuple(t1, o1), std::make_tuple(t2, o2), t1, 1200.),
- o1);
+ std::make_tuple(t1, o1), std::make_tuple(t2, o2),
+ t1 + chrono::nanoseconds(1200), 0.0),
+ o1 + chrono::nanoseconds(60));
EXPECT_EQ(NoncausalTimestampFilter::InterpolateOffsetRemainder(
- std::make_tuple(t1, o1), std::make_tuple(t2, o2), t1, 1200.),
- 60.);
+ std::make_tuple(t1, o1), std::make_tuple(t2, o2),
+ t1 + chrono::nanoseconds(1200), 0.0),
+ 0.0);
+
+ for (int i = -MaxVelocityRatio::den * MaxVelocityRatio::num * 6;
+ i <
+ MaxVelocityRatio::den * MaxVelocityRatio::num * 6 + (t2 - t1).count();
+ ++i) {
+ monotonic_clock::time_point ta_base = t1;
+ const double ta_orig = static_cast<double>(i) / 3.0;
+ double ta = ta_orig;
+
+ NormalizeTimestamps(&ta_base, &ta);
+ CHECK_GE(ta, 0.0);
+ CHECK_LT(ta, 1.0);
+
+ const chrono::nanoseconds expected_offset =
+ NoncausalTimestampFilter::InterpolateOffset(
+ std::make_tuple(t1, o1), std::make_tuple(t2, o2), ta_base);
+
+ EXPECT_EQ(expected_offset, NoncausalTimestampFilter::InterpolateOffset(
+ std::make_tuple(t1, o1),
+ std::make_tuple(t2, o2), ta_base, ta));
+
+ const double expected_double_offset =
+ static_cast<double>(o1.count()) +
+ static_cast<double>(ta_orig) / static_cast<double>((t2 - t1).count()) *
+ (o2 - o1).count();
+
+ EXPECT_NEAR(
+ static_cast<double>(
+ NoncausalTimestampFilter::InterpolateOffset(
+ std::make_tuple(t1, o1), std::make_tuple(t2, o2), ta_base, ta)
+ .count()) +
+ NoncausalTimestampFilter::InterpolateOffsetRemainder(
+ std::make_tuple(t1, o1), std::make_tuple(t2, o2), ta_base, ta),
+ expected_double_offset, 1e-9)
+ << ": i " << i << " t " << ta_base << " " << ta << " t1 " << t1
+ << " o1 " << o1.count() << "ns t2 " << t2 << " o2 " << o2.count()
+ << "ns Non-rounded: " << expected_offset.count() << "ns";
+ }
}
// Tests that all variants of ExtrapolateOffset do reasonable things.
@@ -986,21 +1030,19 @@
const monotonic_clock::time_point t1 = e + chrono::nanoseconds(10000);
const chrono::nanoseconds o1 = chrono::nanoseconds(100);
- //const double o1d = static_cast<double>(o1.count());
const monotonic_clock::time_point t2 = t1 + chrono::nanoseconds(100000);
const chrono::nanoseconds o2 = chrono::nanoseconds(150);
- //const double o2d = static_cast<double>(o2.count());
- EXPECT_EQ(NoncausalTimestampFilter::BoundOffset(
- std::make_tuple(t1, o1), std::make_tuple(t2, o2), t1),
+ EXPECT_EQ(NoncausalTimestampFilter::BoundOffset(std::make_tuple(t1, o1),
+ std::make_tuple(t2, o2), t1),
o1);
EXPECT_EQ(NoncausalTimestampFilter::BoundOffset(
std::make_tuple(t1, o1), std::make_tuple(t2, o2), t1, 0.0),
std::pair(o1, 0.0));
- EXPECT_EQ(NoncausalTimestampFilter::BoundOffset(
- std::make_tuple(t1, o1), std::make_tuple(t2, o2), t2),
+ EXPECT_EQ(NoncausalTimestampFilter::BoundOffset(std::make_tuple(t1, o1),
+ std::make_tuple(t2, o2), t2),
o2);
EXPECT_EQ(NoncausalTimestampFilter::BoundOffset(
std::make_tuple(t1, o1), std::make_tuple(t2, o2), t2, 0.0),
@@ -1013,7 +1055,9 @@
// trust it.
for (int i = -MaxVelocityRatio::den * MaxVelocityRatio::num * 6;
- i < MaxVelocityRatio::den * MaxVelocityRatio::num * 6 + (t2 - t1).count(); ++i) {
+ i <
+ MaxVelocityRatio::den * MaxVelocityRatio::num * 6 + (t2 - t1).count();
+ ++i) {
monotonic_clock::time_point ta_base = t1;
const double ta_orig = static_cast<double>(i) / 3.0;
double ta = ta_orig;
@@ -1238,17 +1282,17 @@
// Confirm the problem statement is reasonable... We've had enough trouble
// here in the past.
- EXPECT_TRUE(
- filter_a.ValidateSolution(&filter_b, Pointer(), t1_a, t1_a + o1_a + chrono::nanoseconds(1)));
- EXPECT_TRUE(
- filter_a.ValidateSolution(&filter_b, Pointer(), t2_a, t2_a + o2_a + chrono::nanoseconds(1)));
+ EXPECT_TRUE(filter_a.ValidateSolution(&filter_b, Pointer(), t1_a,
+ t1_a + o1_a + chrono::nanoseconds(1)));
+ EXPECT_TRUE(filter_a.ValidateSolution(&filter_b, Pointer(), t2_a,
+ t2_a + o2_a + chrono::nanoseconds(1)));
- EXPECT_TRUE(
- filter_b.ValidateSolution(&filter_a, Pointer(), t1_b, t1_b + o1_b + chrono::nanoseconds(1)));
- EXPECT_TRUE(
- filter_b.ValidateSolution(&filter_a, Pointer(), t2_b, t2_b + o2_b + chrono::nanoseconds(1)));
- EXPECT_TRUE(
- filter_b.ValidateSolution(&filter_a, Pointer(), t3_b, t3_b + o3_b + chrono::nanoseconds(1)));
+ EXPECT_TRUE(filter_b.ValidateSolution(&filter_a, Pointer(), t1_b,
+ t1_b + o1_b + chrono::nanoseconds(1)));
+ EXPECT_TRUE(filter_b.ValidateSolution(&filter_a, Pointer(), t2_b,
+ t2_b + o2_b + chrono::nanoseconds(1)));
+ EXPECT_TRUE(filter_b.ValidateSolution(&filter_a, Pointer(), t3_b,
+ t3_b + o3_b + chrono::nanoseconds(1)));
// Before the start
result = filter_a.FindTimestamps(&filter_b, Pointer(),
@@ -1310,15 +1354,12 @@
// doesn't modify the timestamps.
const BootTimestamp t1 = e + chrono::nanoseconds(1000);
const BootDuration o1{0, chrono::nanoseconds(100)};
- const double o1d = static_cast<double>(o1.duration.count());
const BootTimestamp t2 = e + chrono::microseconds(2000);
const BootDuration o2{0, chrono::nanoseconds(150)};
- const double o2d = static_cast<double>(o2.duration.count());
const BootTimestamp t3 = e + chrono::microseconds(3000);
const BootDuration o3{0, chrono::nanoseconds(50)};
- const double o3d = static_cast<double>(o3.duration.count());
const BootTimestamp t4 = e + chrono::microseconds(4000);
@@ -1362,7 +1403,7 @@
filter.Offset(nullptr, Pointer(),
e + (t2.time_since_epoch() + t1.time_since_epoch()) / 2,
0.0, 0),
- std::make_pair(o1, (o2d - o1d) / 2.));
+ std::make_pair(o1 + (o2 - o1) / 2, 0.0));
EXPECT_EQ(filter.Offset(nullptr, Pointer(), t2, 0.0, 0),
std::make_pair(o2, 0.0));
@@ -1371,7 +1412,7 @@
filter.Offset(nullptr, Pointer(),
e + (t2.time_since_epoch() + t3.time_since_epoch()) / 2,
0.0, 0),
- std::make_pair(o2, (o2d + o3d) / 2. - o2d));
+ std::make_pair((o2 + o3) / 2, 0.0));
EXPECT_EQ(filter.Offset(nullptr, Pointer(), t3, 0.0, 0),
std::make_pair(o3, 0.0));