BE-857: atlas: take validated scalars as math operands

libs/@local/graph/atlas/benches/math_kernels.rs

@@ -48,7 +48,7 @@
     clippy::integer_division_remainder_used,
     clippy::significant_drop_tightening,
     reason = "benchmark fixtures compute deterministic floating-point inputs, and Criterion owns \
-              group drops; the crate-level expectations in lib.rs do not extend to bench targets"
+              group drops. The crate-level expectations in lib.rs do not extend to bench targets"
 )]
 
 use core::{hint::black_box, time::Duration};
@@ -327,7 +327,8 @@ fn hardware_counter(event: &str) -> Criterion<darwin_kperf_criterion::HardwareCo
 
     Criterion::default()
         .with_measurement(
-            counter.expect("hardware counters require root on Apple Silicon (run under sudo)"),
+            counter
+                .expect("hardware counter initialization, configuration, and start must succeed"),
         )
         .warm_up_time(Duration::from_millis(500))
         .measurement_time(Duration::from_secs(1))

libs/@local/graph/atlas/src/cli/tui/render/map.rs

@@ -35,7 +35,7 @@ const MAP_MARGIN: Positive = Positive::new(1.04).unwrap();
 ///
 /// The first refresh tick of a collapsed initialization reports rows that all sit together, which
 /// has a centre but no size to scale.
-const MINIMUM_EXTENT: f32 = 1.0;
+const MINIMUM_EXTENT: Positive = Positive::new(1.0).unwrap();
 
 /// Draws the placement itself, the sampled rows as braille dots with the skeleton picked out.
 ///
@@ -52,10 +52,15 @@ pub(super) fn render_map(frame: &mut Frame, area: Rect, placement: &PlacementMap
         .title_bottom(population(placement, area.width).right_aligned());
     let inner = block.inner(area);
 
+    let Some([horizontal, vertical]) = map_bounds(&placement.positions, inner) else {
+        // leave an unrepresentable viewport empty instead of drawing a partial placement.
+        frame.render_widget(block, area);
+        return;
+    };
+
     let (skeleton, interior) = placement.positions.split_at(placement.landmarks);
     let skeleton: Vec<(f64, f64)> = drawable(skeleton).into_iter().collect();
     let interior: Vec<(f64, f64)> = drawable(interior).into_iter().collect();
-    let [horizontal, vertical] = map_bounds(&placement.positions, inner);
 
     let canvas = Canvas::default()
         .block(block)
@@ -96,35 +101,37 @@ fn drawable(positions: &[Vec2]) -> impl IntoIterator<Item = (f64, f64)> {
 /// The map's viewport, which is the placement's own extent squared against the pane's dot grid.
 ///
 /// A braille cell is [`DOTS_ACROSS`] dots wide and [`DOTS_DOWN`] tall over a terminal cell about
-/// twice as tall as it is wide, so a dot is approximately square. Equal data units per dot on both
-/// axes is what keeps the atlas its own shape instead of a version stretched to fill the frame. The
-/// extent grows to the grid's aspect first and then by [`MAP_MARGIN`]. A placement with no extent
-/// of its own grows to [`MINIMUM_EXTENT`] instead, so its rows sit in the middle of a frame rather
-/// than dividing by zero.
-pub(super) fn map_bounds(positions: &[Vec2], inner: Rect) -> [[f64; 2]; 2] {
+/// twice as tall as it is wide, which leaves each dot approximately square. Growing to the grid's
+/// aspect ratio before applying [`MAP_MARGIN`] gives both axes equal data units per dot in exact
+/// arithmetic.
+/// A collapsed placement first grows to [`MINIMUM_EXTENT`] to provide a positive extent.
+///
+/// Returns [`None`] when growth requires a corner beyond the finite `f32` range. With no finite
+/// position, returns a unit viewport around the origin. Corner rounding follows
+/// [`Bounds2::with_aspect_ratio`] and can change the achieved aspect ratio.
+pub(super) fn map_bounds(positions: &[Vec2], inner: Rect) -> Option<[[f64; 2]; 2]> {
     let Some(bounds) = Bounds2::from_points(
         positions
             .iter()
             .copied()
             .filter(|position| position.is_finite()),
     ) else {
-        // Nothing placeable to draw yet.
-        return [[-1.0, 1.0], [-1.0, 1.0]];
+        return Some([[-1.0, 1.0], [-1.0, 1.0]]);
     };
 
     let across = (f32::from(inner.width) * f32::from(DOTS_ACROSS)).max(1.0);
     let down = (f32::from(inner.height) * f32::from(DOTS_DOWN)).max(1.0);
     let aspect = Positive::new(across / down).unwrap_or(Positive::ONE);
 
     let viewport = bounds
-        .with_minimum_extent(MINIMUM_EXTENT)
-        .with_aspect_ratio(aspect)
-        .scaled_about_centre(MAP_MARGIN);
+        .with_minimum_extent(MINIMUM_EXTENT)?
+        .with_aspect_ratio(aspect)?
+        .scaled_about_centre(MAP_MARGIN)?;
 
-    [
+    Some([
         [f64::from(viewport.min().x()), f64::from(viewport.max().x())],
         [f64::from(viewport.min().y()), f64::from(viewport.max().y())],
-    ]
+    ])
 }
 
 /// The map's footer, showing how many rows it is drawing and how many of them are the skeleton.

libs/@local/graph/atlas/src/cli/tui/render/tests.rs

@@ -498,7 +498,8 @@ fn the_map_keeps_the_placement_square() {
         [Vec2::new(-8.0, -1.0), Vec2::new(8.0, 1.0)],
         [Vec2::new(-0.5, -12.0), Vec2::new(0.5, 12.0)],
     ] {
-        let [horizontal, vertical] = map_bounds(&placement, inner);
+        let [horizontal, vertical] =
+            map_bounds(&placement, inner).expect("should represent an ordinary placement viewport");
 
         let across = (horizontal[1] - horizontal[0]) / (f64::from(inner.width) * 2.0);
         let down = (vertical[1] - vertical[0]) / (f64::from(inner.height) * 4.0);
@@ -546,11 +547,13 @@ fn a_placement_with_no_extent_still_has_a_viewport() {
 
     // A collapsed placement and a frame with nothing finite in it
     // both draw a box rather than a degenerate one.
-    let [horizontal, vertical] = map_bounds(&[Vec2::new(3.0, 3.0); 4], inner);
+    let [horizontal, vertical] =
+        map_bounds(&[Vec2::new(3.0, 3.0); 4], inner).expect("should widen the collapsed placement");
     assert!(horizontal[0] < horizontal[1], "{horizontal:?}");
     assert!(vertical[0] < vertical[1], "{vertical:?}");
 
-    let [horizontal, vertical] = map_bounds(&[Vec2::new(f32::NAN, 0.0)], inner);
+    let [horizontal, vertical] =
+        map_bounds(&[Vec2::new(f32::NAN, 0.0)], inner).expect("should preserve the empty viewport");
     assert!(horizontal[0] < horizontal[1], "{horizontal:?}");
     assert!(vertical[0] < vertical[1], "{vertical:?}");
 }

libs/@local/graph/atlas/src/file/generation/tests.rs

@@ -27,8 +27,8 @@ use crate::{
     },
     integrity::{Sha256, Sha256Digest, Update as _},
     math::{
-        AffinityCurve, Bounds2, Vec2, d_non_negative, d_positive, non_negative, open_unit_fraction,
-        unit_fraction,
+        AffinityCurve, Bounds2, Vec2, d_non_negative, d_positive, non_negative, nz,
+        open_unit_fraction, positive, unit_fraction,
     },
     morton::Depth,
     salt::{
@@ -84,8 +84,7 @@ fn config(seed: u64) -> FitConfig {
             maximum_count: NonZero::new(2).expect("the fixture capacity is nonzero"),
             ..
         },
-        curve: AffinityCurve::new(1.577, 0.895)
-            .expect("the fixture parameters are finite and strictly positive"),
+        curve: AffinityCurve::new(positive!(1.577), positive!(0.895)),
         ..
     }
 }

libs/@local/graph/atlas/src/file/salt/tests.rs

@@ -97,8 +97,7 @@ fn config() -> FitConfig {
             maximum_count: NonZero::new(4_096).expect("the fixture capacity is nonzero"),
             ..
         },
-        curve: AffinityCurve::new(1.577, 0.895)
-            .expect("the fixture parameters are finite and strictly positive"),
+        curve: AffinityCurve::new(positive!(1.577), positive!(0.895)),
         placement: placement(),
         policy: PolicyOptions {
             overrides: vec![PolicyOverride {

libs/@local/graph/atlas/src/lib.rs

@@ -142,6 +142,7 @@
 #![allow(
     unused_crate_dependencies,
     unused_features,
+    unused_macros,
     dead_code,
     unreachable_pub,
     unused_imports,

libs/@local/graph/atlas/src/math/affinity/fit.rs

@@ -13,7 +13,9 @@
 use core::num::NonZero;
 
 use super::AffinityCurve;
-use crate::math::{DNonNegative, DPositive, Positive, positive, scalar::narrow_f32};
+use crate::math::{
+    DFinite, DNonNegative, DPositive, Derivation, Positive, d_finite, d_positive, positive,
+};
 
 /// Sample count and distance range for the least-squares target.
 ///
@@ -140,7 +142,7 @@ impl AffinityCurve {
             }
         })?;
 
-        Self::new(narrow_f32(a.get())?, narrow_f32(b.get())?)
+        Some(Self::new(a.narrow()?, b.narrow()?))
     }
 }
 
@@ -193,38 +195,50 @@ impl SampleGrid {
 #[derive(Debug, Copy, Clone)]
 struct NormalEquations {
     /// Sum of squared residuals, the objective the fit minimizes.
-    residual_sum_of_squares: f64,
+    residual_sum_of_squares: DFinite,
     /// The `a`-`a` entry of the normal matrix.
-    j_aa: f64,
+    j_aa: DFinite,
     /// The symmetric off-diagonal entry of the normal matrix.
-    j_ab: f64,
+    j_ab: DFinite,
     /// The `b`-`b` entry of the normal matrix.
-    j_bb: f64,
-    /// The `a` component of the gradient.
-    g_a: f64,
-    /// The `b` component of the gradient.
-    g_b: f64,
+    j_bb: DFinite,
+    /// The `a` component of Jᵀr.
+    g_a: DFinite,
+    /// The `b` component of Jᵀr.
+    g_b: DFinite,
 }
 
-impl NormalEquations {
-    /// The additive identity every accumulation pass starts from.
+/// Unvalidated sums for one objective and Jacobian evaluation.
+struct NormalEquationsDerivation {
+    residual_sum_of_squares: Derivation<DFinite>,
+    j_aa: Derivation<DFinite>,
+    j_ab: Derivation<DFinite>,
+    j_bb: Derivation<DFinite>,
+    g_a: Derivation<DFinite>,
+    g_b: Derivation<DFinite>,
+}
+
+impl NormalEquationsDerivation {
+    /// Empty sums before the first sample.
     const ZERO: Self = Self {
-        residual_sum_of_squares: 0.0,
-        j_aa: 0.0,
-        j_ab: 0.0,
-        j_bb: 0.0,
-        g_a: 0.0,
-        g_b: 0.0,
+        residual_sum_of_squares: Derivation::ZERO,
+        j_aa: Derivation::ZERO,
+        j_ab: Derivation::ZERO,
+        j_bb: Derivation::ZERO,
+        g_a: Derivation::ZERO,
+        g_b: Derivation::ZERO,
     };
 
-    /// Returns whether every accumulated sum is finite.
-    const fn is_finite(self) -> bool {
-        self.residual_sum_of_squares.is_finite()
-            && self.j_aa.is_finite()
-            && self.j_ab.is_finite()
-            && self.j_bb.is_finite()
-            && self.g_a.is_finite()
-            && self.g_b.is_finite()
+    /// Validates the accumulated sums, returning [`None`] if any is non-finite.
+    fn finish(self) -> Option<NormalEquations> {
+        Some(NormalEquations {
+            residual_sum_of_squares: self.residual_sum_of_squares.finish().ok()?,
+            j_aa: self.j_aa.finish().ok()?,
+            j_ab: self.j_ab.finish().ok()?,
+            j_bb: self.j_bb.finish().ok()?,
+            g_a: self.g_a.finish().ok()?,
+            g_b: self.g_b.finish().ok()?,
+        })
     }
 }
 
@@ -334,13 +348,16 @@ fn evaluate(
     a: DPositive,
     b: DPositive,
 ) -> Option<NormalEquations> {
-    let mut sums = NormalEquations::ZERO;
+    let mut sums = NormalEquationsDerivation::ZERO;
+    let exponent = (d_positive!(2.0) * b).finish().ok()?;
 
     for index in 0..grid.samples {
         let distance = grid.distance(index);
-        let power = distance.powf(2.0 * b).get();
-        let denominator = a.get().mul_add(power, 1.0);
-        let residual = 1.0 / denominator - target(distance);
+        let power = distance.powf(exponent.into());
+
+        let denominator = Derivation::from(DNonNegative::from(a)).mul_add(power, DPositive::ONE);
+        let residual = Derivation::from(DFinite::ONE) / denominator - target(distance);
+
         sums.residual_sum_of_squares = residual.mul_add(residual, sums.residual_sum_of_squares);
 
         // the zero-distance sample contributes residual error with zero parameter partials
@@ -349,16 +366,17 @@ fn evaluate(
         };
 
         let denominator_squared = denominator * denominator;
-        let partial_a = -power / denominator_squared;
-        let partial_b = -(2.0 * a * power * distance.ln()) / denominator_squared;
+        let partial_a = (Derivation::from(-DFinite::ONE) * power) / denominator_squared;
+        let partial_b = (d_finite!(-2.0) * a * power * distance.ln()) / denominator_squared;
+
         sums.j_aa = partial_a.mul_add(partial_a, sums.j_aa);
         sums.j_ab = partial_a.mul_add(partial_b, sums.j_ab);
         sums.j_bb = partial_b.mul_add(partial_b, sums.j_bb);
         sums.g_a = partial_a.mul_add(residual, sums.g_a);
         sums.g_b = partial_b.mul_add(residual, sums.g_b);
     }
 
-    sums.is_finite().then_some(sums)
+    sums.finish()
 }
 
 /// Solves the multiplicatively damped 2x2 normal system.
@@ -372,23 +390,56 @@ fn evaluate(
 /// is [`f64::EPSILON`], or when a computed step is non-finite. This numerical floor rejects
 /// near-cancellation, without certifying exact conditioning.
 fn solve_damped(equations: NormalEquations, damping: f64) -> Option<(f64, f64)> {
-    let damped_aa = equations.j_aa * (1.0 + damping);
-    let damped_bb = equations.j_bb * (1.0 + damping);
-    let determinant = damped_aa.mul_add(damped_bb, -(equations.j_ab * equations.j_ab));
-
-    // The damped matrix is positive definite in exact arithmetic; at or
-    // below the floor the closed form divides cancellation noise.
-    if !determinant.is_finite() || determinant <= f64::EPSILON * damped_aa * damped_bb {
+    let damped_aa = Derivation::from(equations.j_aa) * (1.0 + damping);
+    let damped_bb = Derivation::from(equations.j_bb) * (1.0 + damping);
+    let determinant = damped_aa
+        .mul_add(damped_bb, -(equations.j_ab * equations.j_ab))
+        .finish()
+        .ok()?
+        .positive()?;
+
+    // compare determinant cancellation against the product scale of both damped diagonals
+    let determinant_floor = (damped_aa * d_positive!(f64::EPSILON) * damped_bb)
+        .finish()
+        .ok()?;
+    if DFinite::from(determinant) <= determinant_floor {
         return None;
     }
 
-    let step_a = equations
-        .j_ab
+    let step_a = Derivation::from(equations.j_ab)
         .mul_add(equations.g_b, -(damped_bb * equations.g_a))
         / determinant;
-    let step_b = equations
-        .j_ab
+    let step_b = Derivation::from(equations.j_ab)
         .mul_add(equations.g_a, -(damped_aa * equations.g_b))
         / determinant;
-    (step_a.is_finite() && step_b.is_finite()).then_some((step_a, step_b))
+    Some((step_a.finish().ok()?.get(), step_b.finish().ok()?.get()))
+}
+
+#[cfg(test)]
+mod tests {
+    use super::{SampleGrid, evaluate};
+    use crate::math::{DPositive, d_positive};
+
+    #[test]
+    fn evaluation_power_overflow() {
+        // The affinity rounds to zero, but the parameter partials contain ∞/∞.
+        let equations = evaluate(
+            SampleGrid::new(2, d_positive!(1e200)),
+            &|_| 0.0,
+            DPositive::ONE,
+            DPositive::ONE,
+        );
+        assert!(equations.is_none());
+    }
+
+    #[test]
+    fn evaluation_exponent_overflow() {
+        let equations = evaluate(
+            SampleGrid::new(2, DPositive::ONE),
+            &|_| 0.0,
+            DPositive::ONE,
+            d_positive!(f64::MAX),
+        );
+        assert!(equations.is_none());
+    }
 }