mockingbirdnest · Feb 4, 2018 · Feb 4, 2018 · Feb 4, 2018 · Feb 4, 2018
Showing with 44 additions and 13 deletions.

+17 −0 physics/continuous_trajectory.hpp

+27 −13 physics/continuous_trajectory_body.hpp
diff --git a/physics/continuous_trajectory.hpp b/physics/continuous_trajectory.hpp
@@ -1,6 +1,7 @@
 
 #pragma once
 
+#include <atomic>
 #include <experimental/optional>
 #include <utility>
 #include <vector>
@@ -193,6 +194,22 @@ class ContinuousTrajectory : public Trajectory<Frame> {
   // The polynomials are in increasing time order.
   InstantPolynomialPairs polynomials_;
 
+  // Lookups into |polynomials_| are expensive because they entail a binary
+  // search into a vector that grows over time.  In benchmarks, this can be as
+  // costly as the polynomial evaluation itself.  The accesses are not random,
+  // though, they are clustered in time and (slowly) increasing.  To take
+  // advantage of this, we keep track of the index of the last accessed
+  // polynomial and first try to see if the new lookup is for the same
+  // polynomial.  This makes us O(1) instead of O(Log N) most of the time and it
+  // speeds up the lookup by a factor of 7.  This member is atomic because of
+  // multithreading in the ephemeris, and is mutable to maintain the fiction
+  // that evaluation has no side effects.  In the presence of multithreading it
+  // may be that different threads would want to access polynomials at different
+  // indices, but by and large the threads progress in parallel, and benchmarks
+  // show that there is no adverse performance effects.  Any value in the range
+  // of |polynomials_| or 0 is correct.
+  mutable std::atomic_int last_accessed_polynomial_ = 0;
+
   // The time at which this trajectory starts.  Set for a nonempty trajectory.
   std::experimental::optional<Instant> first_time_;
 

diff --git a/physics/continuous_trajectory_body.hpp b/physics/continuous_trajectory_body.hpp
@@ -57,7 +57,7 @@ bool ContinuousTrajectory<Frame>::empty() const {
 
 template<typename Frame>
 double ContinuousTrajectory<Frame>::average_degree() const {
-  if (empty()) {
+  if (polynomials_.empty()) {
     return 0;
   } else {
     double total = 0;
@@ -130,22 +130,24 @@ void ContinuousTrajectory<Frame>::ForgetBefore(Instant const& time) {
   if (polynomials_.empty()) {
     first_time_ = std::experimental::nullopt;
     last_points_.clear();
+    last_accessed_polynomial_ = 0;
   } else {
     first_time_ = time;
+    last_accessed_polynomial_ = polynomials_.size() - 1;
   }
 }
 
 template<typename Frame>
 Instant ContinuousTrajectory<Frame>::t_min() const {
-  if (empty()) {
+  if (polynomials_.empty()) {
     return astronomy::InfiniteFuture;
   }
   return *first_time_;
 }
 
 template<typename Frame>
 Instant ContinuousTrajectory<Frame>::t_max() const {
-  if (empty()) {
+  if (polynomials_.empty()) {
     return astronomy::InfinitePast;
   }
   return polynomials_.crbegin()->t_max;
@@ -452,16 +454,28 @@ template<typename Frame>
 typename ContinuousTrajectory<Frame>::InstantPolynomialPairs::const_iterator
 ContinuousTrajectory<Frame>::FindPolynomialForInstant(
     Instant const& time) const {
-  // Need to use |lower_bound|, not |upper_bound|, because it allows
-  // heterogeneous arguments.  This returns the first polynomial |p| such that
-  // |time <= p.t_max()|.
-  auto const it = std::lower_bound(
-                      polynomials_.begin(), polynomials_.end(), time,
-                      [](InstantPolynomialPair const& left,
-                         Instant const& right) {
-                        return left.t_max < right;
-                      });
-  return it;
+  // This returns the first polynomial |p| such that |time <= p.t_max|.
+  {
+    auto const begin = polynomials_.begin();
+    auto const it = begin + last_accessed_polynomial_;
+    int const last_accessed_polynomial = last_accessed_polynomial_;
+    if (it != polynomials_.end() && time <= it->t_max &&
+        (it == begin || std::prev(it)->t_max < time)) {
+      return it;
+    }
+  }
+  {
+    auto const it =
+        std::lower_bound(polynomials_.begin(),
+                         polynomials_.end(),
+                         time,
+                         [](InstantPolynomialPair const& left,
+                            Instant const& right) {
+                           return left.t_max < right;
+                         });
+    last_accessed_polynomial_ = it - polynomials_.begin();
+    return it;
+  }
 }
 
 }  // namespace internal_continuous_trajectory