[VariationModel] use banker_round when computing deltas

Fixes #1043 and #235

This commit also adds a test from fonttools Tests/varLib/models_test.py that demostrates how our previous approach of late-rounding of deltas can actually lead to an excessive accummulation of rounding errors when the deltas are actually used for interpolation.

fontir/src/variations.rs

@@ -3,7 +3,7 @@ use std::{
     cmp::Ordering,
     collections::{BTreeMap, HashMap, HashSet},
     fmt::{Debug, Display},
-    ops::{Mul, Sub},
+    ops::{Add, Mul, Sub},
 };
 
 use fontdrasil::{
@@ -21,6 +21,44 @@ use write_fonts::{
 
 use crate::error::VariationModelError;
 
+/// Trait for performing banker's rounding on a value.
+///
+/// Banker's rounding rounds half to even. This means that half-way values
+/// are rounded to the nearest even number. For example, 2.5 rounds to 2.0,
+/// 3.5 rounds to 4.0, and -2.5 rounds to -2.0.
+///
+/// This is the same rounding mode as [`f64::round_ties_even`], which was
+/// stabilized in Rust 1.77.0. It also matches Python's buit-in `round` function.
+/// It is usually preferred over alternative approaches because it reduces bias
+/// and errors in calculations.
+///
+/// This trait provides a generic way to apply banker's rounding to different types,
+/// such as `f64` and `kurbo::Vec2`.
+pub trait BankerRound<U, T = Self> {
+    fn banker_round(self) -> U;
+}
+
+impl BankerRound<f64> for f64 {
+    #[inline]
+    fn banker_round(self) -> f64 {
+        self.round_ties_even()
+    }
+}
+
+impl BankerRound<f32> for f32 {
+    #[inline]
+    fn banker_round(self) -> f32 {
+        self.round_ties_even()
+    }
+}
+
+impl BankerRound<kurbo::Vec2> for kurbo::Vec2 {
+    #[inline]
+    fn banker_round(self) -> kurbo::Vec2 {
+        kurbo::Vec2::new(self.x.round_ties_even(), self.y.round_ties_even())
+    }
+}
+
 const ZERO: OrderedFloat<f32> = OrderedFloat(0.0);
 const ONE: OrderedFloat<f32> = OrderedFloat(1.0);
 
@@ -171,7 +209,7 @@ impl VariationModel {
     ) -> Result<Vec<(VariationRegion, Vec<V>)>, DeltaError>
     where
         P: Copy + Default + Sub<P, Output = V>,
-        V: Copy + Mul<f64, Output = V> + Sub<V, Output = V>,
+        V: Copy + Mul<f64, Output = V> + Sub<V, Output = V> + BankerRound<V>,
     {
         if point_seqs.is_empty() {
             return Ok(Vec::new());
@@ -234,7 +272,12 @@ impl VariationModel {
                         })
                         .fold(initial_vector, |acc, (other, other_weight)| {
                             acc - *other * other_weight.into()
-                        }),
+                        })
+                        // deltas will be stored as integers in the VarStore hence must be rounded at
+                        // some point; this is the correct place to round them, instead of at the end,
+                        // otherwise rounding errors can compound especially where master influences
+                        // overlap.
+                        .banker_round(),
                 );
             }
             model_idx_to_result_idx.insert(model_idx, result.len());
@@ -243,6 +286,41 @@ impl VariationModel {
 
         Ok(result)
     }
+
+    /// Convert relative deltas to absolute values at the given location.
+    ///
+    /// Rust version of <https://github.com/fonttools/fonttools/blob/4ad6b0db/Lib/fontTools/varLib/models.py#L514-L545>
+    pub fn interpolate_from_deltas<V>(
+        &self,
+        location: &NormalizedLocation,
+        deltasets: &[(VariationRegion, Vec<V>)],
+    ) -> Vec<V>
+    where
+        V: Copy + Mul<f64, Output = V> + Add<V, Output = V>,
+    {
+        let mut result = None;
+        for (region, deltas) in deltasets.iter() {
+            let scalar = region.scalar_at(location).into_inner() as f64;
+            if scalar == 0.0 {
+                continue;
+            }
+            let contribution = deltas.iter().map(|d| *d * scalar).collect::<Vec<_>>();
+            if result.is_none() {
+                result = Some(contribution);
+            } else {
+                result = Some(
+                    result
+                        .unwrap()
+                        .iter()
+                        .zip(contribution.into_iter())
+                        .map(|(r, c)| *r + c)
+                        .collect(),
+                );
+            }
+        }
+
+        result.unwrap()
+    }
 }
 
 #[derive(Error, Debug)]
@@ -778,6 +856,7 @@ mod tests {
             "ital" => ("Italic", "ital", 0, 0, 1),
             "foo" => ("Foo", "foo ", -1, 0, 1),
             "bar" => ("Bar", "bar ", -1, 0, 1),
+            "axis" => ("Axis", "axis", 0, 0, 1),
             _ => panic!("No definition for {tag}, add it?"),
         };
         let min = UserCoord::new(min as f32);
@@ -1314,14 +1393,123 @@ mod tests {
             (min_wght, vec![5.0]),
         ]);
 
+        let deltas = model.deltas(&point_seqs).unwrap();
+
         assert_eq!(
             vec![
                 (region(&[("wght", 0.0, 0.0, 0.0)]), vec![10.0]),
                 (region(&[("wght", -1.0, -1.0, 0.0)]), vec![-5.0]),
                 (region(&[("wght", 0.0, 1.0, 1.0)]), vec![2.0]),
             ],
-            model.deltas(&point_seqs).unwrap()
+            deltas
         );
+
+        let loc = NormalizedLocation::for_pos(&[("wght", -0.5)]);
+        let expected = vec![7.5];
+        assert_eq!(expected, model.interpolate_from_deltas(&loc, &deltas));
+    }
+
+    #[derive(Debug, Default, Copy, Clone, PartialEq)]
+    struct NoRoundF64(f64);
+
+    impl BankerRound<NoRoundF64> for NoRoundF64 {
+        #[inline]
+        fn banker_round(self) -> NoRoundF64 {
+            self
+        }
+    }
+
+    impl std::ops::Sub for NoRoundF64 {
+        type Output = Self;
+
+        fn sub(self, rhs: Self) -> Self::Output {
+            NoRoundF64(self.0 - rhs.0)
+        }
+    }
+
+    impl std::ops::Mul<f64> for NoRoundF64 {
+        type Output = Self;
+
+        fn mul(self, rhs: f64) -> Self::Output {
+            NoRoundF64(self.0 * rhs)
+        }
+    }
+
+    #[test]
+    fn modeling_error_within_tolerance() {
+        // Compare interpolating from un-rounded float deltas vs rounded deltas,
+        // and check that rounding errors don't accummulate but stay within <= 0.5.
+        // This test was ported from:
+        // https://github.com/fonttools/fonttools/blob/3b9a73ff/Tests/varLib/models_test.py#L167
+        let num_locations = 31;
+        let num_samples = 251;
+        let locations: Vec<_> = (0..num_locations + 1)
+            .map(|i| NormalizedLocation::for_pos(&[("axis", i as f32 / num_locations as f32)]))
+            .collect();
+        let master_values: HashMap<_, _> = locations
+            .iter()
+            .zip((0..num_locations + 1).map(|i| if i == 0 { 0.0 } else { 100.0 }))
+            .map(|(loc, value)| (loc.clone(), vec![value]))
+            .collect();
+        // Same as master_values, but using special f64s that won't get rounded.
+        let master_values_noround: HashMap<_, _> = master_values
+            .iter()
+            .map(|(loc, values)| (loc.clone(), values.iter().map(|v| NoRoundF64(*v)).collect()))
+            .collect();
+
+        let model =
+            VariationModel::new(locations.into_iter().collect(), vec![axis("axis")]).unwrap();
+
+        let mut num_bad_errors = 0;
+        for i in 0..num_samples {
+            let loc = NormalizedLocation::for_pos(&[("axis", i as f32 / num_samples as f32)]);
+
+            // unrounded float deltas
+            let deltas_float: Vec<_> = model
+                .deltas(&master_values_noround)
+                .unwrap()
+                .into_iter()
+                .map(|(region, deltas)| (region, deltas.iter().map(|d| d.0).collect::<Vec<_>>()))
+                .collect();
+            // deltas rounded within the delta computation loop
+            let deltas_round = model.deltas(&master_values).unwrap();
+            // float deltas only rounded at the very end. This is how NOT to round deltas.
+            let deltas_late_round: Vec<_> = deltas_float
+                .iter()
+                .map(|(region, deltas)| {
+                    (
+                        region.clone(),
+                        deltas
+                            .iter()
+                            .map(|d| d.round_ties_even())
+                            .collect::<Vec<_>>(),
+                    )
+                })
+                .collect();
+            // Sanity checks
+            assert_ne!(deltas_float, deltas_round);
+            assert_ne!(deltas_float, deltas_late_round);
+            assert_ne!(deltas_round, deltas_late_round);
+            assert!(deltas_round
+                .iter()
+                .all(|(_, deltas)| { deltas.iter().all(|d| d.fract() == 0.0) }));
+
+            let expected: Vec<f64> = model.interpolate_from_deltas(&loc, &deltas_float);
+            let actual: Vec<f64> = model.interpolate_from_deltas(&loc, &deltas_round);
+
+            let err = (actual[0] - expected[0]).abs();
+            assert!(err <= 0.5, "i: {}, err: {}", i, err);
+
+            // when the influence of many masters overlap a particular location
+            // interpolating from late-rounded deltas may lead to an accummulation
+            // of rounding errors that exceed the max tolerance
+            let bad = model.interpolate_from_deltas(&loc, &deltas_late_round);
+            let err_bad = (bad[0] - expected[0]).abs();
+            if err_bad > 0.5 {
+                num_bad_errors += 1;
+            }
+        }
+        assert!(num_bad_errors > 0);
     }
 
     #[test]