Simplified scalar type to just use num_dual

This commit is contained in:
connor committed 2026-10-02 19:44:00 -04:00
1 parent 6f215788e5
commit e25d654ef1
5 files changed
+31 -87

No files matched your search

+1 -6
View File
@@ -9,20 +9,15 @@ repository = "https://git.rcjohnstone.com/connor/ordinary-diffeq"
license = "MIT"
exclude = ["roadmap/", "hooks/"]
[features]
default = ["autodiff"]
autodiff = ["dep:num-dual"]
[dependencies]
serde = { version = "1.0", features = ["derive"] }
nalgebra = { version = "0.35", features = ["serde-serialize"] }
num-traits = "0.2.19"
num-dual = { version = "0.15", optional = true }
num-dual = "0.15"
roots = "0.0.8"
[dev-dependencies]
approx = "0.5"
num-dual = "0.15"
criterion = "0.7.0"
serde_json = "1.0"
+7 -4
View File
@@ -218,10 +218,13 @@ holds several methods for comparison over the same problem.
## Automatic differentiation
State, time and step size are generic over the `Scalar` trait. `f64` implements it by default,
and so do the forward-mode dual numbers of [`num-dual`](https://crates.io/crates/num-dual) (the
`autodiff` feature, on by default). Solving at a dual scalar propagates exact derivatives alongside
the state.
State, time and step size are generic over the `Scalar` trait, which is a blanket extension of
`DualNum` and `DualStruct` from [`num-dual`](https://crates.io/crates/num-dual): those two supply
the arithmetic, the transcendentals and `re`, which discards derivative information, and `Scalar`
adds only `lift`. Both are implemented for `f64` as well as for every one of that crate's
forward-mode dual numbers, so `Scalar` covers all of them, and a crate depending on this one
alongside another built on `num-dual` can satisfy both with one bound. Solving at a dual scalar
propagates exact derivatives alongside the state.
The only change needed is to the signature of the right-hand side. Here is the pendulum from above,
before and after:
-1
View File
@@ -406,7 +406,6 @@ mod tests {
/// integrator it wraps does, which is only true if nothing erases the
/// scalar. Vern7 is the demanding case: interpolating inside a step calls
/// the right-hand side again, in dual arithmetic, after the solve.
#[cfg(feature = "autodiff")]
#[test]
fn test_method_preserves_autodiff() {
use num_dual::Dual64;
+3 -4
View File
@@ -1,8 +1,7 @@
#![allow(dead_code)]
// The readme is the crate's front page, and including it here means every
// example in it is compiled and run by `cargo test`. Requires `autodiff`
// because one of those examples differentiates.
#![cfg_attr(feature = "autodiff", doc = include_str!("../readme.md"))]
// example in it is compiled and run by `cargo test`.
#![doc = include_str!("../readme.md")]
pub mod callback;
pub mod controller;
@@ -160,7 +159,7 @@ mod tests {
}
/// Tests that the derivatives a dual instantiation produces are the right ones.
#[cfg(all(test, feature = "autodiff"))]
#[cfg(test)]
mod autodiff_tests {
use crate::ode::ProblemFunction;
use crate::prelude::*;
+20 -72
View File
@@ -7,84 +7,33 @@
//! differencing and no extra right-hand-side evaluations beyond the ones the
//! dual arithmetic performs anyway.
//!
//! [`Scalar`] is implemented for `f64` always, and for the dual numbers of the
//! [`num_dual`] crate under the `autodiff` feature (on by default).
//! [`Scalar`] extends the vocabulary [`num_dual`] already defines rather than
//! restating it: [`DualNum`] carries the arithmetic, the transcendentals and
//! `From<f64>`, and [`DualStruct`] carries [`re`](DualStruct::re). Both are
//! implemented for `f64` as well as for every dual type, so one blanket impl
//! covers all of them, and a crate depending on this one alongside another
//! built on `num-dual` can satisfy both with a single bound.
use nalgebra::RealField;
use num_dual::{DualNum, DualStruct};
use num_traits::NumOps;
/// A scalar that an ODE's state, time and step size may take.
pub trait Scalar: RealField + Copy + NumOps<f64, Self> {
/// This value with all derivative information discarded.
fn re(&self) -> f64;
///
/// `re`, which discards derivative information, is inherited from
/// [`DualStruct`]; this trait adds only [`lift`](Scalar::lift).
pub trait Scalar:
RealField + Copy + NumOps<f64, Self> + DualNum<Primitive = f64> + DualStruct<Real = f64>
{
/// A real constant as a value of this type, with zero derivative.
fn lift(x: f64) -> Self;
}
impl Scalar for f64 {
#[inline(always)]
fn re(&self) -> f64 {
*self
}
#[inline(always)]
fn lift(x: f64) -> Self {
x
Self::from(x)
}
}
#[cfg(feature = "autodiff")]
mod dual {
use super::Scalar;
use nalgebra::{Const, RealField};
use num_dual::{DualNum, DualStruct};
/// Implements `Scalar` for one of `num-dual`'s dual number types.
macro_rules! impl_scalar {
($ty:ident) => {
impl<T> Scalar for num_dual::$ty<T>
where
Self: RealField + Copy + DualNum<Primitive = f64> + DualStruct<Real = f64>,
{
#[inline(always)]
fn re(&self) -> f64 {
DualStruct::re(self)
}
#[inline(always)]
fn lift(x: f64) -> Self {
Self::from(x)
}
}
};
($ty:ident, $($n:ident),+) => {
impl<T: nalgebra::Scalar, $(const $n: usize),+> Scalar
for num_dual::$ty<T, $(Const<$n>),+>
where
Self: RealField + Copy + DualNum<Primitive = f64> + DualStruct<Real = f64>,
{
#[inline(always)]
fn re(&self) -> f64 {
DualStruct::re(self)
}
#[inline(always)]
fn lift(x: f64) -> Self {
Self::from(x)
}
}
};
}
impl_scalar!(Dual);
impl_scalar!(Dual2);
impl_scalar!(Dual3);
impl_scalar!(HyperDual);
impl_scalar!(HyperHyperDual);
impl_scalar!(DualVec, N);
impl_scalar!(Dual2Vec, N);
impl_scalar!(HyperDualVec, M, N);
impl<T> Scalar for T where
T: RealField + Copy + NumOps<f64, T> + DualNum<Primitive = f64> + DualStruct<Real = f64>
{
}
#[cfg(test)]
@@ -94,13 +43,12 @@ mod tests {
#[test]
fn test_f64_is_a_scalar() {
assert_eq!(f64::lift(3.5), 3.5);
assert_eq!(Scalar::re(&3.5), 3.5);
assert_eq!(3.5f64.re(), 3.5);
}
#[cfg(feature = "autodiff")]
#[test]
fn test_dual_projects_to_its_real_part() {
use num_dual::{Dual64, DualNum};
use num_dual::Dual64;
// A lifted constant has no derivative; a seeded variable does, and both
// project back to the same real value.
@@ -109,7 +57,7 @@ mod tests {
assert_eq!(c.eps, 0.0);
let x = Dual64::from_re(3.5).derivative();
assert_eq!(Scalar::re(&x), 3.5);
assert_eq!(x.re(), 3.5);
assert_eq!(x.eps, 1.0);
// d/dx x^3 = 3x^2 = 36.75 at x = 3.5