[timekeeper] Final Kalman Filter changes in prep for frequency.

* Receive samples by reference rather than value (since they also need
  to go to the frequency estimator).
* Return the change in estimate from accepting a new sample.
* Use the frequency when building a transform.
* Define a new function to modify the frequency.

Bug: 56868
Change-Id: I48c8e799ec942e7d5c9a4586bbe2045758685eb1
Reviewed-on: https://fuchsia-review.googlesource.com/c/fuchsia/+/514761
Reviewed-by: Satsuki Ueno <satsukiu@google.com>
Commit-Queue: Jody Sankey <jsankey@google.com>
diff --git a/src/sys/time/timekeeper/src/estimator/kalman_filter.rs b/src/sys/time/timekeeper/src/estimator/kalman_filter.rs
index 7f2eb73..400c427 100644
--- a/src/sys/time/timekeeper/src/estimator/kalman_filter.rs
+++ b/src/sys/time/timekeeper/src/estimator/kalman_filter.rs
@@ -3,7 +3,10 @@
 // found in the LICENSE file.
 
 use {
-    crate::{estimator::OSCILLATOR_ERROR_STD_DEV_PPM, time_source::Sample},
+    crate::{
+        estimator::{frequency_to_adjust_ppm, OSCILLATOR_ERROR_STD_DEV_PPM},
+        time_source::Sample,
+    },
     anyhow::{anyhow, Error},
     fuchsia_zircon as zx,
     lazy_static::lazy_static,
@@ -71,12 +74,11 @@
 
 impl KalmanFilter {
     /// Construct a new KalmanFilter initialized to the supplied sample.
-    pub fn new(sample: Sample) -> Self {
-        let Sample { utc, monotonic, std_dev } = sample;
-        let covariance_00 = duration_to_f64(std_dev).powf(2.0).max(MIN_COVARIANCE);
+    pub fn new(Sample { utc, monotonic, std_dev }: &Sample) -> Self {
+        let covariance_00 = duration_to_f64(*std_dev).powf(2.0).max(MIN_COVARIANCE);
         KalmanFilter {
-            reference_utc: utc,
-            monotonic,
+            reference_utc: utc.clone(),
+            monotonic: monotonic.clone(),
             estimate_0: 0f64,
             estimate_1: 1f64,
             covariance_00,
@@ -84,9 +86,9 @@
     }
 
     /// Propagate the estimate forward to the requested monotonic time.
-    fn predict(&mut self, monotonic: zx::Time) {
-        let monotonic_step = duration_to_f64(monotonic - self.monotonic);
-        self.monotonic = monotonic;
+    fn predict(&mut self, monotonic: &zx::Time) {
+        let monotonic_step = duration_to_f64(*monotonic - self.monotonic);
+        self.monotonic = monotonic.clone();
         // Estimated UTC increases by (change in monotonic time) * frequency.
         self.estimate_0 += self.estimate_1 * monotonic_step;
         // Estimated covariance increases as a function of the time step and oscillator error.
@@ -94,9 +96,9 @@
     }
 
     /// Correct the estimate by incorporating measurement data.
-    fn correct(&mut self, utc: zx::Time, std_dev: zx::Duration) {
-        let measurement_variance = duration_to_f64(std_dev).powf(2.0);
-        let measurement_utc_offset = duration_to_f64(utc - self.reference_utc);
+    fn correct(&mut self, utc: &zx::Time, std_dev: &zx::Duration) {
+        let measurement_variance = duration_to_f64(*std_dev).powf(2.0);
+        let measurement_utc_offset = duration_to_f64(*utc - self.reference_utc);
         // Gain is based on the relative variance of the apriori estimate and the new measurement...
         let k_0 = self.covariance_00 / (self.covariance_00 + measurement_variance);
         // ...and determines how much the measurement impacts the apriori estimate...
@@ -105,12 +107,16 @@
         self.covariance_00 = ((1f64 - k_0) * self.covariance_00).max(MIN_COVARIANCE);
     }
 
-    /// Update the estimate to include the supplied sample.
-    pub fn update(&mut self, Sample { utc, monotonic, std_dev }: Sample) -> Result<(), Error> {
+    /// Update the estimate to include the supplied sample, returning the change in estimated UTC
+    /// that this caused.
+    pub fn update(
+        &mut self,
+        Sample { utc, monotonic, std_dev }: &Sample,
+    ) -> Result<zx::Duration, Error> {
         // Ignore any updates that are earlier than the current filter state. Samples from a single
         // time source should arrive in order due to the validation in time_source_manager, but its
         // not impossible that a backwards step occurs during a time source switch.
-        if monotonic < self.monotonic {
+        if *monotonic < self.monotonic {
             return Err(anyhow!(
                 "sample monotonic={} prior to previous monotonic={}",
                 monotonic.into_nanos(),
@@ -120,9 +126,22 @@
 
         // Calculate apriori by moving the estimate forward to the measurement's monotonic time.
         self.predict(monotonic);
+        let apriori_utc = self.utc();
         // Then correct to aposteriori by merging in the measurement.
         self.correct(utc, std_dev);
-        Ok(())
+        let aposteriori_utc = self.utc();
+
+        Ok(aposteriori_utc - apriori_utc)
+    }
+
+    /// Updates the frequency used to estimate utc propagation.
+    #[allow(unused)]
+    pub fn update_frequency(&mut self, frequency: f64) {
+        // These infrequent updates are also a good opportunity to bring the reference utc forward
+        // to the current internal state.
+        self.reference_utc += f64_to_duration(self.estimate_0);
+        self.estimate_0 = 0f64;
+        self.estimate_1 = frequency;
     }
 
     /// Returns a `Transform` describing the estimated synthetic time and error as a function
@@ -131,9 +150,7 @@
         Transform {
             monotonic_offset: self.monotonic.into_nanos(),
             synthetic_offset: self.utc().into_nanos(),
-            // TODO(jsankey): Accommodate an oscillator frequency error when implementing the
-            // frequency correction algorithm.
-            rate_adjust_ppm: 0,
+            rate_adjust_ppm: frequency_to_adjust_ppm(self.estimate_1),
             error_bound_at_offset: ERROR_BOUND_FACTOR as u64 * self.covariance_00.sqrt() as u64,
             error_bound_growth_ppm: ERROR_BOUND_FACTOR * OSCILLATOR_ERROR_STD_DEV_PPM as u32,
         }
@@ -169,7 +186,7 @@
 
     #[fuchsia::test]
     fn initialize() {
-        let filter = KalmanFilter::new(Sample::new(TIME_1 + OFFSET_1, TIME_1, STD_DEV_1));
+        let filter = KalmanFilter::new(&Sample::new(TIME_1 + OFFSET_1, TIME_1, STD_DEV_1));
         let transform = filter.transform();
         assert_eq!(
             transform,
@@ -200,49 +217,126 @@
     fn kalman_filter_performance() {
         // Note: The expected outputs for these test inputs have been validated using the time
         // synchronization simulator we created during algorithm development.
-        let mut filter = KalmanFilter::new(Sample::new(
+        let mut filter = KalmanFilter::new(&Sample::new(
             zx::Time::from_nanos(10001_000000000),
             zx::Time::from_nanos(1_000000000),
             zx::Duration::from_millis(50),
         ));
         assert_eq!(filter.reference_utc, zx::Time::from_nanos(10001_000000000));
         assert_near!(filter.estimate_0, 0f64, 1.0);
+        assert_near!(filter.estimate_1, 1f64, 1e-9);
         assert_near!(filter.covariance_00, 2.5e15, 1.0);
 
-        assert!(filter
-            .update(Sample::new(
-                zx::Time::from_nanos(10101_100000000),
-                zx::Time::from_nanos(101_000000000),
-                zx::Duration::from_millis(200),
-            ))
-            .is_ok());
+        assert_eq!(
+            filter
+                .update(&Sample::new(
+                    zx::Time::from_nanos(10101_100000000),
+                    zx::Time::from_nanos(101_000000000),
+                    zx::Duration::from_millis(200),
+                ))
+                .unwrap(),
+            zx::Duration::from_nanos(100_005887335 - 0 - 100_000000000)
+        );
         assert_near!(filter.estimate_0, 100_005887335.0, 1.0);
+        assert_near!(filter.estimate_1, 1f64, 1e-9);
         assert_near!(filter.covariance_00, 2.3549341505449715e15, 1.0);
 
-        assert!(filter
-            .update(Sample::new(
-                zx::Time::from_nanos(10300_900000000),
-                zx::Time::from_nanos(301_000000000),
-                zx::Duration::from_millis(100),
-            ))
-            .is_ok());
-        assert_near!(filter.estimate_0, 299_985642106.0, 1.0);
+        assert_eq!(
+            filter
+                .update(&Sample::new(
+                    zx::Time::from_nanos(10300_900000000),
+                    zx::Time::from_nanos(301_000000000),
+                    zx::Duration::from_millis(100),
+                ))
+                .unwrap(),
+            zx::Duration::from_nanos(299_985642105 - 100_005887335 - 200_000000000)
+        );
+        assert_near!(filter.estimate_0, 299_985642105.0, 1.0);
+        assert_near!(filter.estimate_1, 1f64, 1e-9);
         assert_near!(filter.covariance_00, 1.9119595120463945e15, 1.0);
     }
 
     #[fuchsia::test]
+    fn frequency_change() {
+        let mut filter = KalmanFilter::new(&Sample::new(
+            zx::Time::from_nanos(10001_000000000),
+            zx::Time::from_nanos(1_000000000),
+            zx::Duration::from_millis(50),
+        ));
+        assert_eq!(
+            filter
+                .update(&Sample::new(
+                    zx::Time::from_nanos(10201_000000000),
+                    zx::Time::from_nanos(201_000000000),
+                    zx::Duration::from_millis(50),
+                ))
+                .unwrap(),
+            zx::Duration::from_nanos(0)
+        );
+        assert_eq!(filter.reference_utc, zx::Time::from_nanos(10001_000000000));
+        assert_near!(filter.estimate_0, 200_000000000.0, 1.0);
+        assert_near!(filter.estimate_1, 1f64, 1e-9);
+        assert_near!(filter.covariance_00, 1252245957276901.8, 1.0);
+        assert_eq!(
+            filter.transform(),
+            Transform {
+                monotonic_offset: 201_000000000,
+                synthetic_offset: 10201_000000000,
+                rate_adjust_ppm: 0,
+                error_bound_at_offset: 70774174,
+                error_bound_growth_ppm: 2 * OSCILLATOR_ERROR_STD_DEV_PPM as u32,
+            }
+        );
+
+        // Updating the frequency should move the internal reference forward to the last sample,
+        // but otherwise doesn't change the offsets and errors reported externally.
+        filter.update_frequency(0.9999);
+        assert_eq!(filter.reference_utc, zx::Time::from_nanos(10201_000000000));
+        assert_near!(filter.estimate_0, 0.0, 1.0);
+        assert_near!(filter.estimate_1, 0.9999, 1e-9);
+        assert_near!(filter.covariance_00, 1252245957276901.8, 1.0);
+        assert_eq!(
+            filter.transform(),
+            Transform {
+                monotonic_offset: 201_000000000,
+                synthetic_offset: 10201_000000000,
+                rate_adjust_ppm: -100,
+                error_bound_at_offset: 70774174,
+                error_bound_growth_ppm: 2 * OSCILLATOR_ERROR_STD_DEV_PPM as u32,
+            }
+        );
+
+        // Even though this new sample has the same monotonic to UTC offset as before, based on the
+        // updated frequency the new UTC is higher than the filter was expecting. It causes an
+        // increase relative to the apriori estimate.
+        assert_eq!(
+            filter
+                .update(&Sample::new(
+                    zx::Time::from_nanos(10301_000000000),
+                    zx::Time::from_nanos(301_000000000),
+                    zx::Duration::from_millis(50),
+                ))
+                .unwrap(),
+            zx::Duration::from_nanos(3341316)
+        );
+        assert_near!(filter.estimate_0, 99993341316.6, 1.0);
+        assert_near!(filter.estimate_1, 0.9999, 1e-9);
+        assert_near!(filter.covariance_00, 835329143746618.4, 1.0);
+    }
+
+    #[fuchsia::test]
     fn covariance_minimum() {
-        let mut filter = KalmanFilter::new(Sample::new(TIME_1 + OFFSET_1, TIME_1, ZERO_DURATION));
+        let mut filter = KalmanFilter::new(&Sample::new(TIME_1 + OFFSET_1, TIME_1, ZERO_DURATION));
         assert_eq!(filter.covariance_00, MIN_COVARIANCE);
-        assert!(filter.update(Sample::new(TIME_2 + OFFSET_2, TIME_2, ZERO_DURATION)).is_ok());
+        assert!(filter.update(&Sample::new(TIME_2 + OFFSET_2, TIME_2, ZERO_DURATION)).is_ok());
         assert_eq!(filter.covariance_00, MIN_COVARIANCE);
     }
 
     #[fuchsia::test]
     fn earlier_monotonic_ignored() {
-        let mut filter = KalmanFilter::new(Sample::new(TIME_2 + OFFSET_1, TIME_2, STD_DEV_1));
+        let mut filter = KalmanFilter::new(&Sample::new(TIME_2 + OFFSET_1, TIME_2, STD_DEV_1));
         assert_near!(filter.estimate_0, 0.0, 1.0);
-        assert!(filter.update(Sample::new(TIME_1 + OFFSET_1, TIME_1, STD_DEV_1)).is_err());
+        assert!(filter.update(&Sample::new(TIME_1 + OFFSET_1, TIME_1, STD_DEV_1)).is_err());
         assert_near!(filter.estimate_0, 0.0, 1.0);
     }
 }
diff --git a/src/sys/time/timekeeper/src/estimator/mod.rs b/src/sys/time/timekeeper/src/estimator/mod.rs
index 8c9d97e..219e1b9 100644
--- a/src/sys/time/timekeeper/src/estimator/mod.rs
+++ b/src/sys/time/timekeeper/src/estimator/mod.rs
@@ -22,6 +22,11 @@
 /// control the growth in error bound and bound the allowed frequency estimates.
 const OSCILLATOR_ERROR_STD_DEV_PPM: u64 = 15;
 
+/// Converts a floating point frequency to a rate adjustment in PPM.
+fn frequency_to_adjust_ppm(frequency: f64) -> i32 {
+    ((frequency - 1.0f64) * 1_000_000f64).round() as i32
+}
+
 /// Maintains an estimate of the relationship between true UTC time and monotonic time on this
 /// device, based on time samples received from one or more time sources.
 #[derive(Debug)]
@@ -37,7 +42,7 @@
 impl<D: Diagnostics> Estimator<D> {
     /// Construct a new estimator initialized to the supplied sample.
     pub fn new(track: Track, sample: Sample, diagnostics: Arc<D>) -> Self {
-        let filter = KalmanFilter::new(sample);
+        let filter = KalmanFilter::new(&sample);
         diagnostics.record(Event::KalmanFilterUpdated {
             track,
             monotonic: filter.monotonic(),
@@ -50,7 +55,7 @@
     /// Update the estimate to include the supplied sample.
     pub fn update(&mut self, sample: Sample) {
         let utc = sample.utc;
-        if let Err(err) = self.filter.update(sample) {
+        if let Err(err) = self.filter.update(&sample) {
             warn!("Rejected update: {}", err);
             return;
         }
@@ -107,6 +112,15 @@
     }
 
     #[fuchsia::test]
+    fn frequency_to_adjust_ppm_test() {
+        assert_eq!(frequency_to_adjust_ppm(0.999), -1000);
+        assert_eq!(frequency_to_adjust_ppm(0.999999), -1);
+        assert_eq!(frequency_to_adjust_ppm(1.0), 0);
+        assert_eq!(frequency_to_adjust_ppm(1.000001), 1);
+        assert_eq!(frequency_to_adjust_ppm(1.001), 1000);
+    }
+
+    #[fuchsia::test]
     fn initialize() {
         let diagnostics = Arc::new(FakeDiagnostics::new());
         let estimator = Estimator::new(