test(ddx-core): gate the finite-difference oracle on a measured noise floor (closes #67)

crates/ddx-core/src/test_utils.rs

@@ -502,6 +502,45 @@ impl Conditioning {
     }
 }
 
+/// How much rounding alone can move a central difference of `f` in `wrt` at
+/// `(x0, y0)` with step `h`: an estimate of the difference quotient's noise
+/// floor, or `None` where `f` does not evaluate.
+///
+/// A finite difference divides `f(v + h) − f(v − h)` by `2h`, so any error in
+/// `f` of size `δ` becomes an error of `δ / h` in the quotient. The error that
+/// matters is not `ε · |f|`: an expression can round an *inner* value and pass
+/// that jitter through a steep outer function. `power(3, power((y + x + 2.5) -
+/// y, 2))` does not depend on `y`, but `(y + x) - y` rounds differently for
+/// each `y`, and the outer `3^(b²)` scales that ulp of jitter by `∂f/∂b ≈ 7e6`,
+/// so the quotient reads `5e-5` where the derivative is exactly `0` (#67). A
+/// Richardson gate does not see it: `fd(h)` and `fd(h/2)` carry noise of the
+/// same size.
+///
+/// So the floor is measured rather than modelled: move `v` by a few ulps, see
+/// how far `f` moves beyond what its slope `slope` accounts for, and divide by
+/// `h`. A point whose floor is not small next to the comparison tolerance
+/// cannot tell a correct derivative from a wrong one.
+pub fn fd_noise_floor(f: &Expr, x0: f64, y0: f64, wrt: Var, h: f64, slope: f64) -> Option<f64> {
+    let v0 = match wrt {
+        Var::X => x0,
+        Var::Y => y0,
+    };
+    let at = |v: f64| match wrt {
+        Var::X => eval(f, v, y0),
+        Var::Y => eval(f, x0, v),
+    };
+    let f0 = at(v0)?;
+    let mut worst: f64 = 0.0;
+    for k in 1..=8 {
+        for sign in [-1.0, 1.0] {
+            let v = v0 + sign * k as f64 * 4.0 * f64::EPSILON * v0.abs().max(1.0);
+            let moved = at(v)? - f0 - slope * (v - v0);
+            worst = worst.max(moved.abs());
+        }
+    }
+    Some(worst / h)
+}
+
 /// Do two floats agree to within `atol + rtol · scale`?
 ///
 /// **Tolerance is relative to the computation scale, not the result.** Two

crates/ddx-core/tests/simulation.rs

@@ -55,10 +55,10 @@ use std::io::Write as _;
 use ddx_core::sqlparser::ast::Expr;
 use ddx_core::sqlparser::dialect::GenericDialect;
 use ddx_core::test_utils::{
-    central_diff, divides_by_noise, eval, gen_adversarial_sql, gen_expr, gen_expr_and_wrt,
-    gen_marker_free_stmt, gen_marker_statement, has_residual_marker, max_intermediate_mag,
-    metamorphic_mismatch, min_domain_margin, parse_expr, run_bounded, seeded, try_parse,
-    try_parse_stmt, Rng, Var,
+    central_diff, divides_by_noise, eval, fd_noise_floor, gen_adversarial_sql, gen_expr,
+    gen_expr_and_wrt, gen_marker_free_stmt, gen_marker_statement, has_residual_marker,
+    max_intermediate_mag, metamorphic_mismatch, min_domain_margin, parse_expr, run_bounded, seeded,
+    try_parse, try_parse_stmt, Rng, Var,
 };
 use ddx_core::{ColRef, Ddx, DiffError};
 
@@ -93,6 +93,8 @@ fn fd_failure(rng: &mut Rng, expr_text: &str, d: &Expr, wrt: Var) -> Option<Stri
                                // Max relative gap between fd(h) and fd(h/2) for the difference to count as
                                // "in its convergent regime" and therefore trustworthy as an oracle.
     const RICHARDSON_TOL: f64 = 1e-4;
+    // How far under the tolerance a point's rounding-noise floor must sit.
+    const NOISE_SAFETY: f64 = 4.0;
     // Above this, some intermediate value is too large for f64 to resolve an
     // O(1) perturbation against — the point is unfit for numeric comparison
     // (total cancellation passes Richardson because *both* fd(h) and fd(h/2)
@@ -155,10 +157,17 @@ fn fd_failure(rng: &mut Rng, expr_text: &str, d: &Expr, wrt: Var) -> Option<Stri
         if (fd_h - fd_h2).abs() > RICHARDSON_TOL * fd_h2.abs().max(1.0) {
             continue;
         }
-        comparable += 1;
         // fd(h/2) is the more accurate estimate at a convergent point.
         let fd = fd_h2;
-        if (fd - dv).abs() > ATOL + RTOL * dv.abs().max(fd.abs()) {
+        let tolerance = ATOL + RTOL * dv.abs().max(fd.abs());
+        // Noise gate: skip points where rounding inside f alone could move the
+        // difference by a fair share of the tolerance (#67).
+        match fd_noise_floor(&f, x0, y0, wrt, H / 2.0, fd) {
+            Some(noise) if NOISE_SAFETY * noise < tolerance => {}
+            _ => continue,
+        }
+        comparable += 1;
+        if (fd - dv).abs() > tolerance {
             disagree += 1;
             if first_bad.is_empty() {
                 first_bad = format!(
@@ -833,6 +842,58 @@ fn differentiation_is_linear_and_obeys_the_product_rule() {
     });
 }
 
+/// The noise gate must be narrow too: it skips points where rounding inside
+/// `f` can move a finite difference by a fair share of the tolerance, and
+/// nothing else (#67).
+#[test]
+fn noise_gate_skips_only_rounding_dominated_points() {
+    let ddx = Ddx::new();
+    // The nightly soak's false positive: `f` does not depend on `y`, but
+    // `(y + x) + 2.5 - y` rounds differently for each `y`, and `3^(b²)` scales
+    // that jitter up past the tolerance. The derivative is 0, and correct.
+    let text = "CAST(power(3, power((((y + x) + abs(2.5)) - y), 2)) AS DOUBLE)";
+    let f = parse_expr(text);
+    let floor = fd_noise_floor(&f, 1.029904, 0.859873, Var::Y, 5e-5, 0.0).unwrap();
+    assert!(
+        floor > 1e-5,
+        "rounding alone moves the quotient by {floor:e}"
+    );
+    let d = ddx.differentiate(&f, &ColRef::bare("y")).unwrap();
+    for seed in 0..32 {
+        let mut rng = seeded(5820328895, seed);
+        assert_eq!(fd_failure(&mut rng, text, &d, Var::Y), None, "seed {seed}");
+    }
+
+    // And the oracle keeps its teeth: a wrong derivative of a well-conditioned
+    // expression is still reported.
+    for (text, wrong) in [("sin(x) * y", "sin(x) * y"), ("exp(x) / y", "exp(x)")] {
+        let mut rng = seeded(1, 0);
+        assert!(
+            fd_failure(&mut rng, text, &parse_expr(wrong), Var::X).is_some(),
+            "d/dx {text} = {wrong} must be caught"
+        );
+    }
+
+    // Well-conditioned expressions keep a floor far under the tolerance, so
+    // the oracle still compares them.
+    for healthy in [
+        "x * y",
+        "sin(x) * exp(y)",
+        "power(x, 3) / y",
+        "ln(x + y) * x",
+    ] {
+        let e = parse_expr(healthy);
+        let floor = fd_noise_floor(&e, 1.029904, 0.859873, Var::X, 5e-5, 0.0).unwrap();
+        let slope = central_diff(&e, 1.029904, 0.859873, Var::X, 5e-5).unwrap();
+        let floor_beyond_slope =
+            fd_noise_floor(&e, 1.029904, 0.859873, Var::X, 5e-5, slope).unwrap();
+        assert!(
+            floor_beyond_slope < 1e-7,
+            "`{healthy}`: floor {floor_beyond_slope:e} (before removing the slope, {floor:e})"
+        );
+    }
+}
+
 // ---------------------------------------------------------------------------
 // Soak test — long-running, #[ignore]-d, driven by env vars (see module docs).
 // ---------------------------------------------------------------------------