Merge branch 'release/0.6.1'

Documentation and housekeeping only; the library behaves exactly as 0.6.0 does.
Three rustdoc warnings in the Jacobian comment are gone, so the formulas there
render as code rather than as attempted links, the tree is rustfmt-clean with the
Butcher tableaux and interpolation coefficients marked to keep their hand
alignment, and the git hooks stay out of the package.

Cargo.toml

@@ -1,13 +1,13 @@
 [package]
 name = "ordinary-diffeq"
-version = "0.6.0"
+version = "0.6.1"
 edition = "2021"
 authors = ["Connor Johnstone"]
 description = "A library for solving differential equations based on the DifferentialEquations.jl julia library."
 readme = "readme.md"
 repository = "https://git.rcjohnstone.com/connor/ordinary-diffeq"
 license = "MIT"
-exclude = ["roadmap/"]
+exclude = ["roadmap/", "hooks/"]
 
 [features]
 default = ["autodiff"]

benches/bs3_vs_dp5.rs

@@ -194,7 +194,11 @@ fn bench_interpolation(c: &mut Criterion) {
             let mut problem = Problem::new(ode, bs3, controller);
             let solution = problem.solve();
             // Interpolate at 100 points
-            black_box((0..100).map(|i| solution.interpolate(i as f64 * 0.05)).collect::<Vec<_>>());
+            black_box(
+                (0..100)
+                    .map(|i| solution.interpolate(i as f64 * 0.05))
+                    .collect::<Vec<_>>(),
+            );
         });
     });
 
@@ -206,7 +210,11 @@ fn bench_interpolation(c: &mut Criterion) {
             let mut problem = Problem::new(ode, dp45, controller);
             let solution = problem.solve();
             // Interpolate at 100 points
-            black_box((0..100).map(|i| solution.interpolate(i as f64 * 0.05)).collect::<Vec<_>>());
+            black_box(
+                (0..100)
+                    .map(|i| solution.interpolate(i as f64 * 0.05))
+                    .collect::<Vec<_>>(),
+            );
         });
     });
 

benches/orbit.rs

@@ -1,7 +1,7 @@
 use criterion::{criterion_group, criterion_main, Criterion};
 
-use ordinary_diffeq::prelude::*;
 use nalgebra::Vector6;
+use ordinary_diffeq::prelude::*;
 
 fn bench_orbit(c: &mut Criterion) {
     let mu = 3.98600441500000e14;

benches/simple_1d.rs

@@ -1,7 +1,7 @@
 use criterion::{criterion_group, criterion_main, Criterion};
 
-use ordinary_diffeq::prelude::*;
 use nalgebra::Vector1;
+use ordinary_diffeq::prelude::*;
 
 fn bench_simple_1d(c: &mut Criterion) {
     type Params = (f64,);

benches/vern7_comparison.rs

@@ -172,7 +172,11 @@ fn bench_vern7_interpolation(c: &mut Criterion) {
             let mut problem = Problem::new(ode, vern7, controller);
             let solution = problem.solve();
             // Interpolate at 100 points - first one computes extra stages
-            black_box((0..100).map(|i| solution.interpolate(i as f64 * 0.05)).collect::<Vec<_>>());
+            black_box(
+                (0..100)
+                    .map(|i| solution.interpolate(i as f64 * 0.05))
+                    .collect::<Vec<_>>(),
+            );
         });
     });
 
@@ -183,7 +187,11 @@ fn bench_vern7_interpolation(c: &mut Criterion) {
             let dp45 = DormandPrince45::new().a_tol(tol).r_tol(tol);
             let mut problem = Problem::new(ode, dp45, controller);
             let solution = problem.solve();
-            black_box((0..100).map(|i| solution.interpolate(i as f64 * 0.05)).collect::<Vec<_>>());
+            black_box(
+                (0..100)
+                    .map(|i| solution.interpolate(i as f64 * 0.05))
+                    .collect::<Vec<_>>(),
+            );
         });
     });
 

hooks/pre-commit

@@ -0,0 +1,11 @@
+#!/bin/sh
+# Fast checks only. On a warm tree these are about a second between them; a hook
+# that takes a minute is a hook that gets --no-verify. The doctests are in
+# pre-push instead, since each one links its own binary and they cost ~15s.
+#
+# Install: git config core.hooksPath hooks
+set -e
+
+cargo fmt --check
+cargo clippy --all-targets -- -D warnings
+cargo test --quiet --lib

hooks/pre-push

@@ -0,0 +1,9 @@
+#!/bin/sh
+# Everything the pre-commit hook is too slow to run. The doctests matter here
+# because the readme is the crate's front page, so its examples are compiled.
+#
+# Install: git config core.hooksPath hooks
+set -e
+
+cargo test --quiet
+cargo doc --quiet --no-deps

roadmap/GETTING_STARTED.md

@@ -127,6 +127,27 @@ This is a solid foundation! The roadmap builds on this.
 
 ## Development Tips
 
+### Checks, and the Hooks That Run Them
+
+The tree is rustfmt-clean, clippy-clean under `-D warnings`, and builds its docs
+without warnings. All three stay that way, and a pre-existing warning is not an
+exemption — fix it or mark it deliberately, as the Butcher tableaux are marked
+`#[rustfmt::skip]` so rustfmt leaves their hand alignment alone.
+
+The hooks in `hooks/` run them. Git does not track where it looks for hooks, so
+each clone needs one command:
+
+```sh
+git config core.hooksPath hooks
+```
+
+`pre-commit` is `cargo fmt --check`, `cargo clippy --all-targets -- -D warnings`
+and `cargo test --lib`, about a second between them on a warm tree. `pre-push`
+adds the full `cargo test` and `cargo doc --no-deps`, which cost around 25
+seconds — the doctests are slow because each one links its own binary, and they
+matter because the readme is the crate's front page and its examples are
+compiled.
+
 ### Testing Strategy
 
 Every feature should have:

src/controller.rs

@@ -76,9 +76,7 @@ impl<T: Scalar> Controller<T> for PIController<T> {
             TryStep::Accepted(prev_step, h)
         } else {
             // Reject the stepsize and propose a smaller one for the current step
-            TryStep::NotYetAccepted(
-                prev_step / self.factor_c1.min(factor_11 / self.safety_factor),
-            )
+            TryStep::NotYetAccepted(prev_step / self.factor_c1.min(factor_11 / self.safety_factor))
         }
     }
 
@@ -124,7 +122,6 @@ impl<T: Scalar> Default for PIController<T> {
     }
 }
 
-
 /// PID (Proportional-Integral-Derivative) step size controller
 ///
 /// The PID controller is an advanced adaptive time-stepping controller that provides
@@ -171,20 +168,20 @@ impl<T: Scalar> Default for PIController<T> {
 #[derive(Debug, Clone, Copy)]
 pub struct PIDController<T: Scalar = f64> {
     // PID Coefficients
-    pub p_term: f64,  // Proportional: reaction to current error
-    pub i_term: f64,  // Integral: reaction to error history
-    pub d_term: f64,  // Derivative: reaction to error rate of change
+    pub p_term: f64, // Proportional: reaction to current error
+    pub i_term: f64, // Integral: reaction to error history
+    pub d_term: f64, // Derivative: reaction to error rate of change
 
     // Constraints
-    pub factor_c1: f64,      // 1/min_factor (maximum step decrease)
-    pub factor_c2: f64,      // 1/max_factor (maximum step increase)
-    pub h_max: f64,          // Maximum allowed step size
-    pub safety_factor: f64,  // Safety factor (typically 0.9)
+    pub factor_c1: f64,     // 1/min_factor (maximum step decrease)
+    pub factor_c2: f64,     // 1/max_factor (maximum step increase)
+    pub h_max: f64,         // Maximum allowed step size
+    pub safety_factor: f64, // Safety factor (typically 0.9)
 
     // Error history for PID control
-    pub err_old: f64,     // ε_{n-1}: previous step error
-    pub err_older: f64,   // ε_{n-2}: error two steps ago
-    pub h_old: f64,       // h_{n-1}: previous step size
+    pub err_old: f64,   // ε_{n-1}: previous step error
+    pub err_older: f64, // ε_{n-2}: error two steps ago
+    pub h_old: f64,     // h_{n-1}: previous step size
 
     // Next step guess
     pub next_step_guess: TryStep<T>,
@@ -199,15 +196,15 @@ impl<T: Scalar> Controller<T> for PIDController<T> {
         let factor = if self.err_old <= 0.0 {
             // First step: use only proportional control
             let factor_11 = err.powf(self.p_term);
-            self.factor_c2.max(
-                self.factor_c1.min(factor_11 / self.safety_factor)
-            )
+            self.factor_c2
+                .max(self.factor_c1.min(factor_11 / self.safety_factor))
         } else if self.err_older <= 0.0 {
             // Second step: use PI control (proportional + integral)
             let factor_11 = err.powf(self.p_term);
             let factor_12 = self.err_old.powf(-self.i_term);
             self.factor_c2.max(
-                self.factor_c1.min(factor_11 * factor_12 / self.safety_factor)
+                self.factor_c1
+                    .min(factor_11 * factor_12 / self.safety_factor),
             )
         } else {
             // Full PID control (proportional + integral + derivative)
@@ -220,7 +217,8 @@ impl<T: Scalar> Controller<T> for PIDController<T> {
                 1.0
             };
             self.factor_c2.max(
-                self.factor_c1.min(factor_11 * factor_12 * factor_13 / self.safety_factor)
+                self.factor_c1
+                    .min(factor_11 * factor_12 * factor_13 / self.safety_factor),
             )
         };
 
@@ -230,7 +228,7 @@ impl<T: Scalar> Controller<T> for PIDController<T> {
 
             // Update error history for next step
             self.err_older = self.err_old;
-            self.err_old = err.max(1.0e-4);  // Prevent very small values
+            self.err_old = err.max(1.0e-4); // Prevent very small values
             self.h_old = h.re();
 
             // Apply maximum step size limit
@@ -297,9 +295,9 @@ impl<T: Scalar> PIDController<T> {
             factor_c2: 1.0 / max_factor,
             h_max: h_max.abs(),
             safety_factor,
-            err_old: 0.0,      // No history initially
-            err_older: 0.0,    // No history initially
-            h_old: 0.0,        // No history initially
+            err_old: 0.0,   // No history initially
+            err_older: 0.0, // No history initially
+            h_old: 0.0,     // No history initially
             next_step_guess: TryStep::NotYetAccepted(T::lift(initial_h)),
         }
     }
@@ -317,14 +315,14 @@ impl<T: Scalar> PIDController<T> {
     pub fn for_order(order: usize) -> Self {
         let k_plus_1 = (order + 1) as f64;
         Self::new(
-            0.49 / k_plus_1,  // p_term: proportional
-            0.34 / k_plus_1,  // i_term: integral
-            0.10 / k_plus_1,  // d_term: derivative
-            10.0,              // max_factor
-            0.2,               // min_factor
-            100000.0,          // h_max
-            0.9,               // safety_factor
-            1e-4,              // initial_h
+            0.49 / k_plus_1, // p_term: proportional
+            0.34 / k_plus_1, // i_term: integral
+            0.10 / k_plus_1, // d_term: derivative
+            10.0,            // max_factor
+            0.2,             // min_factor
+            100000.0,        // h_max
+            0.9,             // safety_factor
+            1e-4,            // initial_h
         )
     }
 }
@@ -342,19 +340,18 @@ impl<T: Scalar> Default for PIDController<T> {
     /// For method-specific tuning, use `PIDController::for_order(order)` instead.
     fn default() -> Self {
         Self::new(
-            0.07,      // p_term (proportional)
-            0.04,      // i_term (integral)
-            0.01,      // d_term (derivative)
-            10.0,      // max_factor
-            0.2,       // min_factor
-            100000.0,  // h_max
-            0.9,       // safety_factor
-            1e-4,      // initial_h
+            0.07,     // p_term (proportional)
+            0.04,     // i_term (integral)
+            0.01,     // d_term (derivative)
+            10.0,     // max_factor
+            0.2,      // min_factor
+            100000.0, // h_max
+            0.9,      // safety_factor
+            1e-4,     // initial_h
         )
     }
 }
 
-
 #[cfg(test)]
 mod tests {
     use super::*;
@@ -406,7 +403,7 @@ mod tests {
         // Default uses conservative PID coefficients
         assert_eq!(controller.p_term, 0.07);
         assert_eq!(controller.i_term, 0.04);
-        assert_eq!(controller.d_term, 0.01);  // True PID with derivative term
+        assert_eq!(controller.d_term, 0.01); // True PID with derivative term
     }
 
     #[test]
@@ -443,7 +440,8 @@ mod tests {
 
         for &err in &errors {
             let mut pi_result = <PIController as Controller<f64>>::determine_step(&mut pi, h, err);
-            let mut pid_result = <PIDController as Controller<f64>>::determine_step(&mut pid, h, err);
+            let mut pid_result =
+                <PIDController as Controller<f64>>::determine_step(&mut pid, h, err);
 
             if pi_result.is_accepted() {
                 pi_steps.push(pi_result.extract());
@@ -555,7 +553,8 @@ mod tests {
 
         for &err in &errors {
             let mut pi_result = <PIController as Controller<f64>>::determine_step(&mut pi, h, err);
-            let mut pid_result = <PIDController as Controller<f64>>::determine_step(&mut pid, h, err);
+            let mut pid_result =
+                <PIDController as Controller<f64>>::determine_step(&mut pid, h, err);
 
             if pi_result.is_accepted() {
                 let pi_h_next = pi_result.reset().unwrap().extract();
@@ -574,11 +573,8 @@ mod tests {
 
         // Compute standard deviation of log step size ratios
         let pi_mean: f64 = pi_ratios.iter().sum::<f64>() / pi_ratios.len() as f64;
-        let pi_variance: f64 = pi_ratios
-            .iter()
-            .map(|r| (r - pi_mean).powi(2))
-            .sum::<f64>()
-            / pi_ratios.len() as f64;
+        let pi_variance: f64 =
+            pi_ratios.iter().map(|r| (r - pi_mean).powi(2)).sum::<f64>() / pi_ratios.len() as f64;
         let pi_std = pi_variance.sqrt();
 
         let pid_mean: f64 = pid_ratios.iter().sum::<f64>() / pid_ratios.len() as f64;

src/integrator/bs3.rs

@@ -64,7 +64,7 @@ pub struct BS3<const D: usize> {
 impl<const D: usize> Default for BS3<D>
 where
     BS3<D>: Integrator<D>,
- {
+{
     fn default() -> Self {
         Self::new()
     }
@@ -120,6 +120,10 @@ impl<'a, const D: usize> BS3Integrator<'a> for BS3<D> {
     // Row 2: [1/2]
     // Row 3: [0, 3/4]
     // Row 4: [2/9, 1/3, 4/9]
+    // The tableau is aligned by hand, so rustfmt leaves it alone.
+    #[rustfmt::skip]
+    // The tableau is aligned by hand, so rustfmt leaves it alone.
+    #[rustfmt::skip]
     const A: &'a [f64] = &[
         1.0 / 2.0,           // a[1,0]
         0.0,                 // a[2,0]
@@ -130,6 +134,8 @@ impl<'a, const D: usize> BS3Integrator<'a> for BS3<D> {
     ];
 
     // Solution weights (3rd order)
+    #[rustfmt::skip]
+    #[rustfmt::skip]
     const B: &'a [f64] = &[
         2.0 / 9.0,           // b[0]
         1.0 / 3.0,           // b[1]
@@ -138,6 +144,8 @@ impl<'a, const D: usize> BS3Integrator<'a> for BS3<D> {
     ];
 
     // Error estimate weights (difference between 3rd and 2nd order)
+    #[rustfmt::skip]
+    #[rustfmt::skip]
     const B_ERROR: &'a [f64] = &[
         2.0 / 9.0 - 7.0 / 24.0,      // b[0] - b*[0]
         1.0 / 3.0 - 1.0 / 4.0,       // b[1] - b*[1]
@@ -146,6 +154,8 @@ impl<'a, const D: usize> BS3Integrator<'a> for BS3<D> {
     ];
 
     // Stage times
+    #[rustfmt::skip]
+    #[rustfmt::skip]
     const C: &'a [f64] = &[
         0.0,                 // c[0]
         1.0 / 2.0,           // c[1]
@@ -191,7 +201,9 @@ where
 
         // Stage 4: k4 = f(t + c[3]*h, y + h*(a[3,0]*k1 + a[3,1]*k2 + a[3,2]*k3))
         let y4 = ode.y
-            + (k[0] * T::lift(Self::A[3]) + k[1] * T::lift(Self::A[4]) + k[2] * T::lift(Self::A[5]))
+            + (k[0] * T::lift(Self::A[3])
+                + k[1] * T::lift(Self::A[4])
+                + k[2] * T::lift(Self::A[5]))
                 * h;
         k[3] = (ode.f)(ode.t + h * Self::C[3], y4, &ode.params);
 
@@ -225,10 +237,10 @@ where
         // BS3 uses standard cubic Hermite interpolation with derivatives at endpoints
         // Store: y0, y1, f0=k[0], f1=k[3] (FSAL)
         let dense_coeffs = vec![
-            ode.y,           // y0 at start of step
-            next_y,          // y1 at end of step
-            k[0],            // f(t0, y0) - derivative at start
-            k[3],            // f(t1, y1) - derivative at end (FSAL)
+            ode.y,  // y0 at start of step
+            next_y, // y1 at end of step
+            k[0],   // f(t0, y0) - derivative at start
+            k[3],   // f(t1, y1) - derivative at end (FSAL)
         ];
 
         StepResult {
@@ -304,7 +316,7 @@ mod tests {
         let ode = ODE::new(&derivative, 0.0, 1.0, y0, ());
 
         let bs3 = BS3::new();
-        let h = 0.001;  // Smaller step size for tighter tolerances
+        let h = 0.001; // Smaller step size for tighter tolerances
 
         let __step = bs3.step(&ode, h);
         let (y_next, err, dense) = (__step.y, __step.err, __step.dense);
@@ -334,7 +346,7 @@ mod tests {
         let ode = ODE::new(&derivative, 0.0, 1.0, y0, ());
 
         let bs3 = BS3::new();
-        let h = 0.001;  // Smaller step size
+        let h = 0.001; // Smaller step size
 
         let __step = bs3.step(&ode, h);
         let (_y_next, _err, dense) = (__step.y, __step.err, __step.dense);

src/integrator/dormand_prince.rs

@@ -212,7 +212,11 @@ where
             for (j, item) in k.iter().enumerate().take(i) {
                 y_term += item * T::lift(Self::A[(i * (i - 1)) / 2 + j]);
             }
-            k[i] = (ode.f)(ode.t + T::lift(Self::C[i]) * h, ode.y + y_term * h, &ode.params);
+            k[i] = (ode.f)(
+                ode.t + T::lift(Self::C[i]) * h,
+                ode.y + y_term * h,
+                &ode.params,
+            );
 
             // Use that and bis to calculate the y and error terms
             next_y += k[i] * h * T::lift(Self::B[i]);

src/integrator/rosenbrock.rs

@@ -28,10 +28,11 @@ impl Default for JacobianUpdateStrategy {
 
 /// Compute the Jacobian matrix ∂f/∂y using forward finite differences
 ///
-/// For a system y' = f(t, y), computes the D×D Jacobian matrix J where J[i,j] = ∂f_i/∂y_j
+/// For a system `y' = f(t, y)`, computes the D×D Jacobian matrix `J` where
+/// `J[i,j] = ∂f_i/∂y_j`
 ///
-/// Uses forward differences: J[i,j] ≈ (f_i(y + ε*e_j) - f_i(y)) / ε
-/// where ε = √(machine_epsilon) * max(|y[j]|, 1.0)
+/// Uses forward differences: `J[i,j] ≈ (f_i(y + ε*e_j) - f_i(y)) / ε`
+/// where `ε = √(machine_epsilon) * max(|y[j]|, 1.0)`
 pub fn compute_jacobian<T: Scalar, const D: usize, P>(
     f: ProblemFunction<'_, T, D, P>,
     t: T,
@@ -443,12 +444,21 @@ mod tests {
         // So log(error) = 2*log(h) + const
         // Slope should be approximately 2
         for i in 0..errors.len() - 1 {
-            let rate =
-                (errors[i].log10() - errors[i + 1].log10()) / (step_sizes[i].log10() - step_sizes[i + 1].log10());
-            println!("Convergence rate between h={} and h={}: {}", step_sizes[i], step_sizes[i+1], rate);
+            let rate = (errors[i].log10() - errors[i + 1].log10())
+                / (step_sizes[i].log10() - step_sizes[i + 1].log10());
+            println!(
+                "Convergence rate between h={} and h={}: {}",
+                step_sizes[i],
+                step_sizes[i + 1],
+                rate
+            );
 
             // Should be close to 2 for a 2nd order method
-            assert!(rate > 1.8 && rate < 2.2, "Convergence rate {} is not close to 2", rate);
+            assert!(
+                rate > 1.8 && rate < 2.2,
+                "Convergence rate {} is not close to 2",
+                rate
+            );
         }
     }
 
@@ -466,7 +476,7 @@ mod tests {
         let analytical = |t: f64| 0.5 * (1.01 * t).exp();
 
         // Test convergence with Julia's step sizes: (1/2)^(6:-1:3) = [1/64, 1/32, 1/16, 1/8]
-        let step_sizes = vec![1.0/64.0, 1.0/32.0, 1.0/16.0, 1.0/8.0];
+        let step_sizes = vec![1.0 / 64.0, 1.0 / 32.0, 1.0 / 16.0, 1.0 / 8.0];
         let mut errors = Vec::new();
 
         for &h in &step_sizes {
@@ -504,8 +514,11 @@ mod tests {
         println!("Orders per step refinement: {:?}", orders);
 
         // Julia tests: @test sim.𝒪est[:final]≈2 atol=0.2
-        assert!((avg_order - 2.0).abs() < 0.2,
-                "Convergence order {:.2} not within 0.2 of expected order 2", avg_order);
+        assert!(
+            (avg_order - 2.0).abs() < 0.2,
+            "Convergence order {:.2} not within 0.2 of expected order 2",
+            avg_order
+        );
     }
 
     #[test]
@@ -528,11 +541,16 @@ mod tests {
         let mut problem = Problem::new(ode, rb23, controller);
         let solution = problem.solve();
 
-        println!("Adaptive solve completed in {} steps", solution.states.len());
+        println!(
+            "Adaptive solve completed in {} steps",
+            solution.states.len()
+        );
 
-        assert!(solution.states.len() < 20,
-                "Adaptive solve should complete in less than 20 steps, got {}",
-                solution.states.len());
+        assert!(
+            solution.states.len() < 20,
+            "Adaptive solve should complete in less than 20 steps, got {}",
+            solution.states.len()
+        );
 
         // Verify final value is accurate
         let analytical = 0.5 * (1.01_f64 * 1.0).exp();

src/integrator/vern7.rs

@@ -80,7 +80,7 @@ pub struct Vern7<const D: usize> {
 impl<const D: usize> Default for Vern7<D>
 where
     Vern7<D>: Integrator<D>,
- {
+{
     fn default() -> Self {
         Self::new()
     }
@@ -141,6 +141,8 @@ impl<'a, const D: usize> Vern7Integrator<'a> for Vern7<D> {
     // Stage 8: [a81, 0, a83, a84, a85, a86, a87]
     // Stage 9: [a91, 0, a93, a94, a95, a96, a97, a98]
     // Stage 10: [a101, 0, a103, a104, a105, a106, a107, 0, 0]
+    // The tableau is aligned by hand, so rustfmt leaves it alone.
+    #[rustfmt::skip]
     const A: &'a [f64] = &[
         // Stage 2
         0.005,
@@ -177,6 +179,7 @@ impl<'a, const D: usize> Vern7Integrator<'a> for Vern7<D> {
     ];
 
     // Error estimate weights (difference between 7th and 6th order: b - b*)
+    #[rustfmt::skip]
     const B_ERROR: &'a [f64] = &[
         0.002547011879931045,   // b1 - b*1
         0.0,                    // b2 - b*2
@@ -191,6 +194,7 @@ impl<'a, const D: usize> Vern7Integrator<'a> for Vern7<D> {
     ];
 
     // Time nodes (c coefficients)
+    #[rustfmt::skip]
     const C: &'a [f64] = &[
         0.0,                  // c1
         0.005,                // c2
@@ -210,6 +214,7 @@ impl<'a, const D: usize> Vern7Integrator<'a> for Vern7<D> {
 
 impl<'a, const D: usize> Vern7ExtraStages<'a> for Vern7<D> {
     // Time nodes for extra stages
+    #[rustfmt::skip]
     const C_EXTRA: &'a [f64] = &[
         1.0,    // c11
         0.29,   // c12
@@ -375,6 +380,8 @@ where
         }
     }
 
+    // One interpolation coefficient per line, likewise.
+    #[rustfmt::skip]
     fn interpolate<T: Scalar>(
         &self,
         t_start: T,
@@ -744,7 +751,7 @@ mod tests {
         }
 
         let y0 = Vector1::new(1.0);
-        let t_end: f64 = 1.0;  // Longer interval to get larger errors
+        let t_end: f64 = 1.0; // Longer interval to get larger errors
         let exact = t_end.exp();
 
         let step_sizes: [f64; 3] = [0.2, 0.1, 0.05];
@@ -771,16 +778,8 @@ mod tests {
         let ratio2 = errors[1] / errors[2];
 
         // Allow some tolerance (expect ratio between 64 and 256)
-        assert!(
-            ratio1 > 64.0 && ratio1 < 256.0,
-            "First ratio: {}",
-            ratio1
-        );
-        assert!(
-            ratio2 > 64.0 && ratio2 < 256.0,
-            "Second ratio: {}",
-            ratio2
-        );
+        assert!(ratio1 > 64.0 && ratio1 < 256.0, "First ratio: {}", ratio1);
+        assert!(ratio2 > 64.0 && ratio2 < 256.0, "Second ratio: {}", ratio2);
     }
 
     #[test]
@@ -828,7 +827,7 @@ mod tests {
             Vector2::new(y[1], -y[0])
         }
 
-        let y0 = Vector2::new(1.0, 0.0);  // Start at maximum displacement, zero velocity
+        let y0 = Vector2::new(1.0, 0.0); // Start at maximum displacement, zero velocity
 
         // Period of harmonic oscillator is 2π
         let period = 2.0 * std::f64::consts::PI;
@@ -854,7 +853,10 @@ mod tests {
 
         println!("Initial energy: {}", energy_initial);
         println!("Final energy: {}", energy_final);
-        println!("Energy drift after {} periods: {:.2e}", num_periods, energy_drift);
+        println!(
+            "Energy drift after {} periods: {:.2e}",
+            num_periods, energy_drift
+        );
         println!("Number of steps: {}", solution.times.len());
 
         // Energy should be conserved to high precision (< 1e-7 relative error over 1000 periods)

src/lib.rs

@@ -19,8 +19,8 @@ pub mod prelude {
     pub use super::integrator::method::{Method, MethodName};
     pub use super::integrator::rosenbrock::Rosenbrock23;
     pub use super::integrator::vern7::Vern7;
-    pub use super::ode::ODE;
     pub use super::integrator::{Integrator, StepResult, Tolerances};
+    pub use super::ode::ODE;
     pub use super::problem::{Problem, Solution};
     pub use super::scalar::Scalar;
 }
@@ -376,7 +376,13 @@ mod autodiff_tests {
         }
 
         let y0 = Dual64::from_re(1.0).derivative();
-        let ode = ODE::new(&derivative, Dual64::lift(0.0), Dual64::lift(2.0), Vector1::new(y0), ());
+        let ode = ODE::new(
+            &derivative,
+            Dual64::lift(0.0),
+            Dual64::lift(2.0),
+            Vector1::new(y0),
+            (),
+        );
         let vern7 = Vern7::new().a_tol(1e-12).r_tol(1e-12);
         let mut problem = Problem::new(ode, vern7, PIController::default());
         let solution = problem.solve();

src/problem.rs

@@ -91,8 +91,7 @@ where
         // crossing starts from rather than at — see `EVENT_PROBE_FRACTION`.
         let mut last_fired: Vec<Option<f64>> = vec![None; self.callbacks.len()];
         while !self.ode.finished && self.ode.t.re() < self.ode.t_end.re() {
-            if (self.ode.t + self.controller.next_step_guess().extract()).re()
-                > self.ode.t_end.re()
+            if (self.ode.t + self.controller.next_step_guess().extract()).re() > self.ode.t_end.re()
             {
                 // If the next step would go past the end, then just set it to the end
                 self.controller
@@ -149,7 +148,9 @@ where
                 // The event condition is a real-valued question, so it is
                 // asked in the reals — see this module's header.
                 let before = if search_start.re() == 0.0 {
-                    callback.event(self.ode.t, &self.ode.y, &self.ode.params).re()
+                    callback
+                        .event(self.ode.t, &self.ode.y, &self.ode.params)
+                        .re()
                 } else {
                     let probe = self.integrator.step(&self.ode, search_start).y;
                     callback
@@ -171,7 +172,9 @@ where
                     let f = |test_h: f64| {
                         let h = T::lift(test_h);
                         let test_y = self.integrator.step(&self.ode, h).y;
-                        callback.event(self.ode.t + h, &test_y, &self.ode.params).re()
+                        callback
+                            .event(self.ode.t + h, &test_y, &self.ode.params)
+                            .re()
                     };
                     Some(
                         find_root_brent(search_start.re(), curr_step.re(), &f, &mut convergency)
@@ -366,6 +369,7 @@ mod tests {
     use crate::callback::{stop, FnCallback};
     use crate::controller::PIController;
     use crate::integrator::dormand_prince::DormandPrince45;
+    use crate::integrator::vern7::Vern7;
     use crate::integrator::Tolerances;
     use approx::assert_relative_eq;
     use nalgebra::{Matrix3, Vector1, Vector2, Vector3};
@@ -440,8 +444,7 @@ mod tests {
             effect: &|ode: &mut ODE<'_, f64, 1, ()>| ode.y = Vector1::new(100.0),
         };
 
-        let mut problem =
-            Problem::new(ode, dp45, PIController::default()).with_callback(&mut jump);
+        let mut problem = Problem::new(ode, dp45, PIController::default()).with_callback(&mut jump);
         let solution = problem.solve();
 
         // Five units of integration before the jump, five after.
@@ -464,8 +467,7 @@ mod tests {
             effect: &|ode: &mut ODE<'_, f64, 1, ()>| ode.y = Vector1::new(100.0),
         };
 
-        let mut problem =
-            Problem::new(ode, dp45, PIController::default()).with_callback(&mut jump);
+        let mut problem = Problem::new(ode, dp45, PIController::default()).with_callback(&mut jump);
         let solution = problem.solve();
 
         let at_crossing: Vec<usize> = (0..solution.times.len())
@@ -650,8 +652,7 @@ mod tests {
             effect: &|ode: &mut ODE<'_, f64, 1, ()>| ode.y = -ode.y,
         };
 
-        let mut problem =
-            Problem::new(ode, dp45, PIController::default()).with_callback(&mut flip);
+        let mut problem = Problem::new(ode, dp45, PIController::default()).with_callback(&mut flip);
         let solution = problem.solve();
 
         // y rises to 1 at t = 1 and flips to -1, which the effect leaves on the
@@ -665,7 +666,11 @@ mod tests {
         assert_relative_eq!(jumps[0], 1.0, max_relative = 1e-9);
         assert_relative_eq!(jumps[1], 3.0, max_relative = 1e-9);
         assert_relative_eq!(solution.times.last().unwrap(), &3.0, max_relative = 1e-12);
-        assert_relative_eq!(solution.states.last().unwrap()[0], -1.0, max_relative = 1e-8);
+        assert_relative_eq!(
+            solution.states.last().unwrap()[0],
+            -1.0,
+            max_relative = 1e-8
+        );
     }
 
     /// The controller's error history is discarded when an effect changes the
@@ -708,8 +713,8 @@ mod tests {
                 effect,
             };
             let ode = ODE::new(&unit_rate, 0.0, 10.0, Vector1::new(0.0), ());
-            let mut problem = Problem::new(ode, DormandPrince45::new(), spy)
-                .with_callback(&mut callback);
+            let mut problem =
+                Problem::new(ode, DormandPrince45::new(), spy).with_callback(&mut callback);
             problem.solve();
             resets.get()
         }
@@ -739,7 +744,7 @@ mod tests {
         }
 
         let ode = ODE::new(&free_fall, 0.0, 2.0, Vector2::new(1.0, 0.0), G);
-        let vern7 = crate::integrator::vern7::Vern7::new().a_tol(1e-12).r_tol(1e-12);
+        let vern7 = Vern7::new().a_tol(1e-12).r_tol(1e-12);
         let mut bounce = FnCallback {
             event: &|_t: f64, y: &Vector2<f64>, _p: &f64| y[0],
             effect: &|ode: &mut ODE<'_, f64, 2, f64>| {
@@ -803,7 +808,7 @@ mod tests {
         let omega = Vector3::new(0.0, 0.0, 1.0);
         let t_end = 5.0; // past half a turn, short of a full one
         let ode = ODE::new(&mrp_rate, 0.0, t_end, Vector3::zeros(), omega);
-        let vern7 = crate::integrator::vern7::Vern7::new().a_tol(1e-12).r_tol(1e-12);
+        let vern7 = Vern7::new().a_tol(1e-12).r_tol(1e-12);
 
         let mut shadow = FnCallback {
             event: &|_t: f64, s: &Vector3<f64>, _p: &Vector3<f64>| 1.0 - s.norm(),
@@ -830,7 +835,11 @@ mod tests {
         // And the attitude at the end is the one the body is actually in.
         let expected = -Vector3::new(0.0, 0.0, 1.0 / (t_end / 4.0).tan());
         for i in 0..3 {
-            assert_relative_eq!(solution.states.last().unwrap()[i], expected[i], epsilon = 1e-9);
+            assert_relative_eq!(
+                solution.states.last().unwrap()[i],
+                expected[i],
+                epsilon = 1e-9
+            );
         }
     }
 

src/scalar.rs

@@ -115,5 +115,4 @@ mod tests {
         // d/dx x^3 = 3x^2 = 36.75 at x = 3.5
         assert_eq!(x.powi(3).eps, 3.0 * 3.5 * 3.5);
     }
-
 }