connor/nrlmsise00
1
0 Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Simplified the scalar to use actual num_dual methods

Browse Source
This commit is contained in:
connor committed 2026-10-02 19:40:37 -04:00
1 parent f35faf45bb
commit 169bc7a6bd
9 files changed
+242 -305
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

Cargo.lock
Generated
-1
View File
@@ -329,7 +329,6 @@ version = "0.1.0"
dependencies = [
"approx",
"criterion",
"nalgebra",
"num-dual",
"num-traits",
"serde_json",
Cargo.toml
+8 -9
View File
@@ -25,9 +25,7 @@ all-features = true
members = ["nrlmsis2"]
[features]
default = ["autodiff", "drivers"]
# Dual-number scalars, so density can be differentiated.
autodiff = ["dep:num-dual", "dep:nalgebra"]
default = ["drivers"]
# The compiled-in solar and geomagnetic record, and date lookups into it.
drivers = []
# Trades the reference's exact bits for speed: the angle-sum identity for the expansions'
@@ -36,13 +34,14 @@ drivers = []
fast = []
[dependencies]
nalgebra = { version = "0.35", optional = true }
num-dual = { version = "0.15", optional = true }
num-dual = { version = "0.15", default-features = false }
num-traits = { version = "0.2", default-features = false }
[dev-dependencies]
approx = "0.5"
criterion = "0.7"
# With default features, for the benchmark: it measures the vector dual types, which are
# the ones behind num-dual's nalgebra feature.
num-dual = "0.15"
# arbitrary_precision keeps each number's own digits, so a fixture value is parsed by
# str::parse rather than by serde_json's own float parser, which reads some f64s back one
@@ -53,10 +52,10 @@ serde_json = { version = "1.0", features = ["arbitrary_precision"] }
[[bench]]
name = "density"
harness = false
# Measures the dual-number scalars and the driver lookup, so it needs both.
required-features = ["autodiff", "drivers"]
# Measures the driver lookup as well as the model.
required-features = ["drivers"]
[[example]]
name = "density"
# Reads the compiled-in driver record and seeds a dual number, so it needs both.
required-features = ["drivers", "autodiff"]
# Reads the compiled-in driver record.
required-features = ["drivers"]
readme.md
+45 -8
View File
@@ -12,7 +12,8 @@ implementations exist and both are `f64` only.
use nrlmsise00::{Input, MagneticActivity, Model};
let model = Model::new();
let output = model.gtd7d(&Input::consistent(
// The scalar is named, because a call built only from literals has nothing else to fix it.
let output = model.gtd7d::<f64>(&Input::consistent(
172, // day of year
29000.0, // seconds of universal time
400.0, // altitude, km
@@ -32,8 +33,6 @@ assert!((output.temperature - 1098.2).abs() < 0.1);
And the derivative, through the model rather than around it:
```rust
# #[cfg(feature = "autodiff")]
# fn main() {
use num_dual::{Dual64, DualNum};
use nrlmsise00::{Input, MagneticActivity, Model, Scalar};
@@ -53,15 +52,49 @@ let output = model.gtd7d(&Input::consistent(
// d(rho)/dz, in kg/m^3 per km, and the scale height it implies.
let scale_height = -output.total_mass.re / output.total_mass.eps;
assert!((scale_height - 55.8).abs() < 0.1);
# }
# #[cfg(not(feature = "autodiff"))]
# fn main() {}
```
`cargo run --example density` prints a vertical profile, the same altitude across two
solar cycles, and the altitude derivative a drag calculation needs, checked against a
central difference.
## The scalar type is `num-dual`'s
`Scalar` is a marker over [`num_dual::DualNum`] with one method of its own, `lift`, and it
is blanket-implemented:
```rust,ignore
pub trait Scalar: DualNum<Primitive = f64> + DualStruct<Real = f64> + Copy {
fn lift(x: f64) -> Self { Self::from(x) }
}
impl<T: DualNum<Primitive = f64> + DualStruct<Real = f64> + Copy> Scalar for T {}
```
So every type that implements `num-dual`'s trait at `f64` precision is a scalar here —
`f64`, each of its dual numbers, nestings of them, and anything a caller writes. There is
nothing to implement and nothing to enumerate, which is the point: an earlier version
declared `exp`, `ln`, `powf` and the rest itself and implemented them for a fixed list of
types, so a caller whose scalar was a dual type the list had missed was stopped by the
orphan rule with no way out, and a caller generic over another crate's scalar trait of the
same shape was adapting between two names for one thing.
Two things follow for a caller:
- **Generic code imports only `Scalar`.** The other methods are its supertraits', and a
type parameter bounded by `Scalar` carries them. `DualNum`, `DualStruct` and `Signed` are
re-exported from the `scalar` module for code that calls them on a concrete dual number.
- **A call built only from float literals has to name its scalar**, as
`model.gtd7d::<f64>(..)`. Inference will not settle a literal on `f64` when the blanket
bound admits every dual number as well. Code that gets its position from somewhere typed
never meets this.
`num-dual` is therefore not optional, and this crate takes it with `default-features =
false`: that drops `nalgebra`, which only the vector dual types need and which this crate
never names. A caller who wants `DualSVec64` enables it from their own manifest, and the
blanket implementation covers those types without this crate knowing about them.
## `gtd7` or `gtd7d`
They differ in one number. `gtd7`'s total mass density sums helium, atomic oxygen,
@@ -94,7 +127,9 @@ about the model rather than anyone's code.
below says what it costs.
At a dual type the value agrees to about 1e-14 rather than exactly, because `num-dual`
computes a division as a reciprocal and a multiply. See the `scalar` module docs.
computes a division as a reciprocal and a multiply. At `f64` every method is a thin forward
to the `std` function of the same name, which is why that path is the exact one. See the
`scalar` module docs.
## Where it deliberately differs from the reference
@@ -252,7 +287,9 @@ use nrlmsise00::{drivers, Model};
let model = Model::new();
let when = drivers::mjd(2015, 3, 15, 8, 3, 20.0); // or hifitime's to_mjd_utc_days()
let output = model.gtd7d_at(when, 400.0, 60.0, -70.0).expect("inside the record");
let output = model
.gtd7d_at::<f64>(when, 400.0, 60.0, -70.0)
.expect("inside the record");
assert!((output.total_mass - 2.3e-12).abs() < 0.5e-12);
# }
src/expansion.rs
+13 -14
View File
@@ -548,18 +548,17 @@ pub(crate) fn globe7<T: Scalar>(
* (plg[0][1] * p[119] * swc[5] * cd14 + 1.0)
* ((plg[0][1] * p[68] + plg[0][3] * p[69] + plg[0][5] * p[70])
* shifted(work.universal, phase[P71], SR, input.seconds_of_day, p[71]));
t[11] = t[11]
+ (plg[2][3] * p[76] + plg[2][5] * p[77] + plg[2][7] * p[78])
* swc[11]
* mixed(
work.universal,
work.east2,
phase[P79],
input.seconds_of_day,
p[79],
input.longitude_deg,
)
* (dfa * p[137] * swc[1] + 1.0);
t[11] += (plg[2][3] * p[76] + plg[2][5] * p[77] + plg[2][7] * p[78])
* swc[11]
* mixed(
work.universal,
work.east2,
phase[P79],
input.seconds_of_day,
p[79],
input.longitude_deg,
)
* (dfa * p[137] * swc[1] + 1.0);
}
// Mixed magnetic activity, universal time and longitude.
@@ -621,7 +620,7 @@ pub(crate) fn globe7<T: Scalar>(
let mut tinf = T::lift(p[30]);
for (index, term) in t.iter().enumerate() {
let switch = if index == 8 { sw9 } else { sw[index + 1] };
tinf = tinf + *term * switch.abs();
tinf += *term * switch.abs();
}
tinf
}
@@ -735,7 +734,7 @@ pub(crate) fn glob7s<T: Scalar>(
let mut tt = T::lift(0.0);
for (index, term) in t.iter().enumerate() {
let switch = if index == 8 { sw9 } else { sw[index + 1] };
tt = tt + *term * switch.abs();
tt += *term * switch.abs();
}
tt
}
src/helpers.rs
+10 -11
View File
@@ -225,15 +225,14 @@ pub(crate) fn splini<T: Scalar>(xa: &[T], ya: &[T], y2a: &[T], x: T) -> T {
let b = (xx - xa[klo]) / h;
let a2 = a * a;
let b2 = b * b;
yi = yi
+ ((T::lift(1.0) - a2) * ya[klo] / 2.0
+ b2 * ya[khi] / 2.0
+ ((-(a2 * a2 + 1.0) / 4.0 + a2 / 2.0) * y2a[klo]
+ (b2 * b2 / 4.0 - b2 / 2.0) * y2a[khi])
* h
* h
/ 6.0)
* h;
yi += ((T::lift(1.0) - a2) * ya[klo] / 2.0
+ b2 * ya[khi] / 2.0
+ ((-(a2 * a2 + 1.0) / 4.0 + a2 / 2.0) * y2a[klo]
+ (b2 * b2 / 4.0 - b2 / 2.0) * y2a[khi])
* h
* h
/ 6.0)
* h;
klo += 1;
khi += 1;
}
@@ -271,8 +270,8 @@ mod tests {
/// coefficient puts there.
#[test]
fn test_gravity_varies_with_latitude_as_the_coefficient_says() {
let (equator, _) = glatf(0.0);
let (pole, _) = glatf(90.0);
let (equator, _): (f64, f64) = glatf(0.0);
let (pole, _): (f64, f64) = glatf(90.0);
assert!(pole > equator);
let fraction = (pole - equator) / equator;
// 1 - 0.0026373 * cos(0) against 1 - 0.0026373 * cos(pi), so twice the coefficient.
src/model.rs
+23 -24
View File
@@ -43,7 +43,7 @@ const NUCLEON_MASS: f64 = 1.66e-24;
/// use nrlmsise00::{Input, MagneticActivity, Model};
///
/// let model = Model::new();
/// let output = model.gtd7d(&Input::consistent(
/// let output = model.gtd7d::<f64>(&Input::consistent(
/// 172, // day of year
/// 29000.0, // seconds of universal time
/// 400.0, // km
@@ -138,7 +138,7 @@ impl Model {
let mut output = self.evaluate(input, &mut work);
output.total_mass = mass_density(&output, true);
if self.switches.units == Units::Si {
output.total_mass = output.total_mass / 1000.0;
output.total_mass /= 1000.0;
}
output
}
@@ -363,7 +363,7 @@ impl Model {
// the output units, so only the thousand is applied here.
output.total_mass = mass_density(&output, false);
if metric {
output.total_mass = output.total_mass / 1000.0;
output.total_mass /= 1000.0;
}
// The temperature at the altitude, from the same profile with no mass.
@@ -530,7 +530,7 @@ impl Model {
let rl = (b28 * PDM[0][1] / b04).ln();
let zc04 = PDM[0][4] * PDL[1][0];
let hc04 = PDM[0][5] * PDL[1][1];
output.helium = output.helium * ccor(z, rl, hc04, zc04);
output.helium *= ccor(z, rl, hc04, zc04);
}
// Atomic oxygen.
@@ -563,12 +563,12 @@ impl Model {
let hc16 = PDM[1][5] * PDL[1][3];
let zc16 = PDM[1][4] * PDL[1][2];
let hc216 = PDM[1][5] * PDL[1][4];
output.atomic_oxygen = output.atomic_oxygen * ccor2(z, rl, hc16, zc16, hc216);
output.atomic_oxygen *= ccor2(z, rl, hc16, zc16, hc216);
// Chemistry correction.
let hcc16 = PDM[1][7] * PDL[1][13];
let zcc16 = PDM[1][6] * PDL[1][12];
let rc16 = PDM[1][3] * PDL[1][14];
output.atomic_oxygen = output.atomic_oxygen * ccor(z, T::lift(rc16), hcc16, zcc16);
output.atomic_oxygen *= ccor(z, T::lift(rc16), hcc16, zcc16);
}
// Molecular oxygen.
@@ -600,7 +600,7 @@ impl Model {
let rl = (b28 * PDM[3][1] / b32).ln();
let hc32 = PDM[3][5] * PDL[1][7];
let zc32 = PDM[3][4] * PDL[1][6];
output.molecular_oxygen = output.molecular_oxygen * ccor(z, rl, hc32, zc32);
output.molecular_oxygen *= ccor(z, rl, hc32, zc32);
}
// Correction for a general departure from diffusive equilibrium above the
// lower boundary, which applies at every altitude.
@@ -609,8 +609,7 @@ impl Model {
let zcc32 = PDM[3][6] * PDL[1][21];
let rc32 = ((input.f107_average - 150.0) * (sw[1] * PDL[0][23]) + 1.0)
* (PDM[3][3] * PDL[1][23]);
output.molecular_oxygen =
output.molecular_oxygen * ccor2(z, rc32, hcc32, zcc32, hcc232);
output.molecular_oxygen *= ccor2(z, rc32, hcc32, zcc32, hcc232);
}
// Argon.
@@ -641,7 +640,7 @@ impl Model {
let rl = (b28 * PDM[4][1] / b40).ln();
let hc40 = PDM[4][5] * PDL[1][9];
let zc40 = PDM[4][4] * PDL[1][8];
output.argon = output.argon * ccor(z, rl, hc40, zc40);
output.argon *= ccor(z, rl, hc40, zc40);
}
// Hydrogen.
@@ -673,12 +672,12 @@ impl Model {
let rl = (b28 * PDM[5][1] * PDL[1][17].abs() / b01).ln();
let hc01 = PDM[5][5] * PDL[1][11];
let zc01 = PDM[5][4] * PDL[1][10];
output.hydrogen = output.hydrogen * ccor(z, rl, hc01, zc01);
output.hydrogen *= ccor(z, rl, hc01, zc01);
// Chemistry correction.
let hcc01 = PDM[5][7] * PDL[1][19];
let zcc01 = PDM[5][6] * PDL[1][18];
let rc01 = PDM[5][3] * PDL[1][20];
output.hydrogen = output.hydrogen * ccor(z, T::lift(rc01), hcc01, zcc01);
output.hydrogen *= ccor(z, T::lift(rc01), hcc01, zcc01);
}
// Atomic nitrogen.
@@ -709,12 +708,12 @@ impl Model {
let rl = (b28 * PDM[6][1] * PDL[0][2].abs() / b14).ln();
let hc14 = PDM[6][5] * PDL[0][1];
let zc14 = PDM[6][4] * PDL[0][0];
output.atomic_nitrogen = output.atomic_nitrogen * ccor(z, rl, hc14, zc14);
output.atomic_nitrogen *= ccor(z, rl, hc14, zc14);
// Chemistry correction.
let hcc14 = PDM[6][7] * PDL[0][4];
let zcc14 = PDM[6][6] * PDL[0][3];
let rc14 = PDM[6][3] * PDL[0][5];
output.atomic_nitrogen = output.atomic_nitrogen * ccor(z, T::lift(rc14), hcc14, zcc14);
output.atomic_nitrogen *= ccor(z, T::lift(rc14), hcc14, zcc14);
}
// Anomalous oxygen: hot atomic oxygen, on its own barometric profile rather than
@@ -793,14 +792,14 @@ fn mass_density<T: Scalar>(output: &Output<T>, anomalous: bool) -> T {
/// mass density is divided by a thousand. Both operations, in that order: g/cm^3 scaled by
/// 1e6 and then by 1e-3 is kg/m^3, and doing it in one step would not round the same way.
fn to_si<T: Scalar>(output: &mut Output<T>) {
output.helium = output.helium * 1.0e6;
output.atomic_oxygen = output.atomic_oxygen * 1.0e6;
output.molecular_nitrogen = output.molecular_nitrogen * 1.0e6;
output.molecular_oxygen = output.molecular_oxygen * 1.0e6;
output.argon = output.argon * 1.0e6;
output.total_mass = output.total_mass * 1.0e6;
output.hydrogen = output.hydrogen * 1.0e6;
output.atomic_nitrogen = output.atomic_nitrogen * 1.0e6;
output.anomalous_oxygen = output.anomalous_oxygen * 1.0e6;
output.total_mass = output.total_mass / 1000.0;
output.helium *= 1.0e6;
output.atomic_oxygen *= 1.0e6;
output.molecular_nitrogen *= 1.0e6;
output.molecular_oxygen *= 1.0e6;
output.argon *= 1.0e6;
output.total_mass *= 1.0e6;
output.hydrogen *= 1.0e6;
output.atomic_nitrogen *= 1.0e6;
output.anomalous_oxygen *= 1.0e6;
output.total_mass /= 1000.0;
}
src/profile.rs
+1 -1
View File
@@ -106,7 +106,7 @@ pub(crate) fn densu<T: Scalar>(
if tt.re() <= 0.0 {
expl = T::lift(50.0);
}
densu_temp = dlb * (tlb / tt).powt(gamma + (1.0 + alpha)) * expl;
densu_temp = dlb * (tlb / tt).powd(gamma + (1.0 + alpha)) * expl;
if alt.re() >= za {
return densu_temp;
}
src/scalar.rs
+138 -233
View File
@@ -5,27 +5,90 @@
//! exists: density falls by a factor of e about every 50 km, so a drag partial that drops
//! its dependence on position is badly wrong rather than slightly wrong.
//!
//! [`Scalar`] is implemented for `f64` always, and for the dual numbers of the
//! [`num_dual`] crate under the `autodiff` feature, which is on by default.
//! # It is `num-dual`'s trait, not one of ours
//!
//! [`Scalar`] declares one method of its own and inherits the rest from
//! [`num_dual::DualNum`], which `num-dual` implements for `f64` and `f32` as well as for
//! its dual types. It is then blanket-implemented, so **every** type satisfying that trait
//! at `f64` precision is a scalar here, including nestings and types this crate has never
//! heard of.
//!
//! That is deliberate. An earlier version declared `exp`, `ln`, `powf` and the rest itself
//! and implemented them for a list of types. Two things were wrong with that. The list
//! could not be complete, and nobody else could extend it: a caller whose scalar was a
//! dual type the list had missed was stopped by the orphan rule, with no way to implement
//! this crate's trait for it. And declaring the same methods a second time made this trait
//! incompatible with every other crate's trait over the same types — `ordinary-diffeq`
//! declares one of the same shape — so a caller generic over both was adapting between
//! two names for one thing. Taking the bound from `num-dual` leaves one name.
//!
//! A caller generic over another crate's scalar trait still writes both bounds, since
//! neither trait implies the other, but the second one is now a formality that any dual
//! number satisfies rather than a wall that some do not.
//!
//! # What has to be in scope
//!
//! Nothing, to write generic code. The methods belong to [`Scalar`]'s supertraits, and a
//! type parameter bounded by `Scalar` carries them without the supertraits being imported:
//!
//! ```
//! use nrlmsise00::Scalar;
//!
//! fn scale_height<T: Scalar>(density: T, slope: T) -> T {
//! // ln from DualNum, abs from Signed, re from DualStruct, lift from Scalar, and only
//! // Scalar is imported.
//! let _ = (density.ln(), density.abs(), density.re());
//! -density / slope + T::lift(0.0)
//! }
//!
//! assert!((scale_height(1.0_f64, -0.02) - 50.0).abs() < 1e-12);
//! ```
//!
//! Calling the same methods on a *concrete* dual number does need the trait that declares
//! them, and the three are re-exported here for that: [`DualNum`] for the arithmetic,
//! [`DualStruct`] for [`DualStruct::re`], [`Signed`] for [`Signed::abs`].
//!
//! # A literal-only call has to name its scalar
//!
//! The one wrinkle of the bound being a blanket one. Inference will not settle a float
//! literal on `f64` when the only thing constraining it is `T: Scalar`, because the
//! blanket implementation admits every dual number too, so a call built entirely from
//! literals is ambiguous and says so at the first method call on its result. Name the
//! scalar once and the whole call follows:
//!
//! ```
//! # use nrlmsise00::{Input, MagneticActivity, Model};
//! let model = Model::new();
//! let input = Input::consistent(
//! 172, 29000.0, 400.0, 60.0, -70.0, 150.0, 150.0,
//! MagneticActivity::Daily { ap: 4.0 },
//! );
//! let output = model.gtd7d::<f64>(&input);
//! assert!(output.total_mass > 0.0);
//! ```
//!
//! Code that gets its position and time from somewhere with a type never meets this.
//!
//! # Writing against it
//!
//! Three rules, each of which the model's own code obeys throughout:
//!
//! - **Magnitude comparisons branch on [`Scalar::re`].** A dual number carries a
//! derivative that an ordering cannot account for, and every comparison in the model —
//! the altitude bands, the exponent clamps, the guards in `ccor` and `dnet` — is asking
//! about a real quantity, where the real part is the right thing to test.
//! - **Magnitude comparisons branch on [`DualStruct::re`].** `num-dual` does order its
//! dual numbers, and it does so on the real part alone, so `x > y` would agree with
//! `x.re() > y.re()`. The model writes the second anyway: every comparison in it — the
//! altitude bands, the exponent clamps, the guards in `ccor` and `dnet` — is a question
//! about a real quantity, and saying so at each site is worth the few characters.
//! - **Coefficients stay `f64`.** `T op f64` is available, so lifting three thousand
//! coefficients to `T` would only cost memory and cache.
//! - **`f64` on the left of an operator does not exist.** `T * 2.0` works; `2.0 * T`
//! cannot, because of the orphan rule. Write the scalar on the left, or
//! [`Scalar::lift`] the constant.
//! [`Scalar::lift`] the constant. Where the reference puts a constant on the left of a
//! *non-commutative* operator, lifting is the only way to keep its arithmetic.
//!
//! # A dual evaluation's value is not bit for bit the `f64` one
//!
//! It agrees to about 1e-14 relative. The difference is `num-dual`'s arithmetic rather
//! than anything this crate does, and two places account for it:
//! It agrees to about 1e-14. The difference is `num-dual`'s arithmetic rather than
//! anything this crate does, and two places account for it:
//!
//! - **Division goes through the reciprocal.** `Dual / Dual` is computed as
//! `self.re * other.re.recip()`, two roundings where a division is one, so a dual
@@ -39,245 +102,87 @@
//! Both were measured rather than assumed. Correcting them here was tried and dropped: a
//! real part can be fixed exactly, by adding back a difference that is itself exactly
//! representable, but the same cannot be done for division without replacing every `/` in
//! the model, so the result was a slower model that was still not exact. The `f64` path is
//! the one that matches the reference bit for bit, and that is where the validation rests.
//!
//! # A consumer generic over its own scalar trait
//!
//! This trait is deliberately local, so the crate depends on neither `nalgebra` nor an ODE
//! solver. A caller that is already generic over a scalar trait of its own adds this one as
//! a second bound:
//!
//! ```ignore
//! fn drag<T: my_crate::Scalar + nrlmsise00::Scalar>(...) { ... }
//! ```
//!
//! Both are satisfied by `f64` and by `num-dual`'s types, so nothing else changes.
//! the model, so the result was a slower model that was still not exact. At `f64` every
//! method is a thin forward to the `std` function of the same name, so that path is the
//! one that matches the reference bit for bit, and it is where the validation rests.
use core::fmt::Debug;
use core::ops::{Add, Div, Mul, Neg, Sub};
pub use num_dual::{DualNum, DualStruct};
pub use num_traits::Signed;
/// A scalar the model can be evaluated in.
///
/// The methods are the ones NRLMSISE-00 actually computes with, and no more.
pub trait Scalar:
Copy
+ Debug
+ Add<Output = Self>
+ Sub<Output = Self>
+ Mul<Output = Self>
+ Div<Output = Self>
+ Neg<Output = Self>
+ Add<f64, Output = Self>
+ Sub<f64, Output = Self>
+ Mul<f64, Output = Self>
+ Div<f64, Output = Self>
{
/// This value with all derivative information discarded.
///
/// Every magnitude comparison in the model goes through here.
fn re(&self) -> f64;
/// Satisfied by `f64`, by `num-dual`'s dual numbers over it, and by anything else that
/// implements [`DualNum`] at `f64` precision and is [`Copy`]. There is nothing to
/// implement: the blanket implementation below covers every such type.
///
/// The methods the model computes with are [`DualNum`]'s, [`DualStruct::re`] and
/// [`Signed::abs`]; see the module documentation for what to import.
pub trait Scalar: DualNum<Primitive = f64> + DualStruct<Real = f64> + Copy {
/// A real constant as a value of this type, with zero derivative.
fn lift(x: f64) -> Self;
/// `e^self`.
fn exp(self) -> Self;
/// The natural logarithm.
fn ln(self) -> Self;
/// The base-10 logarithm.
fn log10(self) -> Self;
/// The square root.
fn sqrt(self) -> Self;
/// The sine, in radians.
fn sin(self) -> Self;
/// The cosine, in radians.
fn cos(self) -> Self;
/// The magnitude.
///
/// The reference writes `sqrt(x*x)` for this. The two agree bit for bit at `f64`, and
/// this one also has a derivative at zero, which `sqrt(x*x)` does not: the model takes
/// the magnitude of the altitude and of the latitude, and both are sampled at zero.
fn abs(self) -> Self;
/// `self` raised to an integer power.
fn powi(self, n: i32) -> Self;
/// `self` raised to a real power.
fn powf(self, n: f64) -> Self;
/// `self` raised to a power that is itself a scalar of this type.
///
/// At `f64` this is the platform `pow`, which is what the reference calls.
fn powt(self, n: Self) -> Self;
}
impl Scalar for f64 {
#[inline(always)]
fn re(&self) -> f64 {
*self
}
#[inline(always)]
/// `num-dual` spells this `Self::from`, which is a less obvious thing to read in the
/// middle of an expansion than `T::lift`, and which the model writes a hundred and
/// twenty times.
fn lift(x: f64) -> Self {
x
}
#[inline(always)]
fn exp(self) -> Self {
f64::exp(self)
}
#[inline(always)]
fn ln(self) -> Self {
f64::ln(self)
}
#[inline(always)]
fn log10(self) -> Self {
f64::log10(self)
}
#[inline(always)]
fn sqrt(self) -> Self {
f64::sqrt(self)
}
#[inline(always)]
fn sin(self) -> Self {
f64::sin(self)
}
#[inline(always)]
fn cos(self) -> Self {
f64::cos(self)
}
#[inline(always)]
fn abs(self) -> Self {
f64::abs(self)
}
#[inline(always)]
fn powi(self, n: i32) -> Self {
f64::powi(self, n)
}
#[inline(always)]
fn powf(self, n: f64) -> Self {
f64::powf(self, n)
}
#[inline(always)]
fn powt(self, n: Self) -> Self {
f64::powf(self, n)
Self::from(x)
}
}
#[cfg(feature = "autodiff")]
mod dual {
use super::Scalar;
use num_dual::{DualNum, DualStruct};
impl<T: DualNum<Primitive = f64> + DualStruct<Real = f64> + Copy> Scalar for T {}
/// Implements `Scalar` for one of `num-dual`'s dual number types.
///
/// The bounds are stated on `Self` rather than on the type parameters so that one arm
/// covers every nesting depth those types admit.
macro_rules! impl_scalar {
($ty:ident) => {
impl<T> Scalar for num_dual::$ty<T>
where
Self: Copy + DualNum<Primitive = f64> + DualStruct<Real = f64>,
{
impl_scalar_body!();
}
};
($ty:ident, $($n:ident),+) => {
impl<T: nalgebra::Scalar $(, const $n: usize)+> Scalar
for num_dual::$ty<T $(, nalgebra::Const<$n>)+>
where
Self: Copy + DualNum<Primitive = f64> + DualStruct<Real = f64>,
{
impl_scalar_body!();
}
};
#[cfg(test)]
mod tests {
use super::*;
/// The blanket implementation reaches the primitive, the dual types, and a nesting of
/// them that this crate never names — which is the point of it being blanket.
#[test]
fn test_the_blanket_implementation_reaches_what_it_should() {
fn assert_scalar<T: Scalar>(x: T) -> f64 {
// One method from each of the four traits, to show they are all available.
(x.exp().ln().abs() + T::lift(1.0)).re()
}
assert_eq!(assert_scalar(2.0_f64), 3.0);
assert_eq!(assert_scalar(num_dual::Dual64::from(2.0)), 3.0);
assert_eq!(assert_scalar(num_dual::Dual2_64::from(2.0)), 3.0);
assert_eq!(assert_scalar(num_dual::HyperDual64::from(2.0)), 3.0);
// A dual of a dual, which the old list of implementations did not include.
assert_eq!(
assert_scalar(num_dual::Dual::<num_dual::Dual64>::from(2.0)),
3.0
);
}
/// The method bodies, which are the same for every one of those types.
macro_rules! impl_scalar_body {
() => {
#[inline(always)]
fn re(&self) -> f64 {
DualStruct::re(self)
}
/// `lift` is a constant: the value with no derivative, whatever the order.
#[test]
fn test_lift_carries_no_derivative() {
use num_dual::Dual64;
#[inline(always)]
fn lift(x: f64) -> Self {
Self::from(x)
}
let constant = Dual64::lift(3.5);
assert_eq!(constant.re, 3.5);
assert_eq!(constant.eps, 0.0);
#[inline(always)]
fn exp(self) -> Self {
DualNum::exp(&self)
}
#[inline(always)]
fn ln(self) -> Self {
DualNum::ln(&self)
}
#[inline(always)]
fn log10(self) -> Self {
DualNum::log10(&self)
}
#[inline(always)]
fn sqrt(self) -> Self {
DualNum::sqrt(&self)
}
#[inline(always)]
fn sin(self) -> Self {
DualNum::sin(&self)
}
#[inline(always)]
fn cos(self) -> Self {
DualNum::cos(&self)
}
#[inline(always)]
fn abs(self) -> Self {
num_traits::Signed::abs(&self)
}
#[inline(always)]
fn powi(self, n: i32) -> Self {
DualNum::powi(&self, n)
}
#[inline(always)]
fn powf(self, n: f64) -> Self {
DualNum::powf(&self, n)
}
#[inline(always)]
fn powt(self, n: Self) -> Self {
DualNum::powd(&self, n)
}
};
let seeded = Dual64::from(3.5).derivative();
assert_eq!(seeded.re(), 3.5);
assert_eq!(seeded.eps, 1.0);
}
impl_scalar!(Dual);
impl_scalar!(Dual2);
impl_scalar!(Dual3);
impl_scalar!(HyperDual);
impl_scalar!(DualVec, N);
impl_scalar!(Dual2Vec, N);
/// Ordering a dual number compares its real part, which is what the model's branches
/// want, so the two forms agree.
#[test]
fn test_ordering_a_dual_compares_the_real_part() {
use num_dual::Dual64;
let small = Dual64::from(1.0).derivative();
let large = Dual64::from(2.0);
assert!(small < large);
assert!(small.re() < large.re());
// Equal real parts compare equal even though the derivatives differ, where a
// lexicographic ordering would have made one of them the smaller.
let other = Dual64::from(1.0);
assert_ne!(small.eps, other.eps);
assert_eq!(small.partial_cmp(&other), Some(core::cmp::Ordering::Equal));
}
}
tests/derivatives.rs
+4 -4
View File
@@ -1,5 +1,3 @@
#![cfg(feature = "autodiff")]
//! That the model differentiates.
//!
//! This is the one thing a port buys that depending on an existing `f64` implementation
@@ -116,7 +114,9 @@ fn test_the_model_jumps_at_the_species_threshold_altitudes() {
// join between the Bates profile above and the spline below, and an eighth at
// 428.475 km, too small to pin here, is where atomic nitrogen's turbopause blend
// saturates and `dnet` switches to the mixed density with a step of about exp(-10).
let thresholds = [
// Annotated, because `Scalar` is blanket-implemented over `num-dual`'s trait and a
// bare float literal flowing into the generic model no longer falls back to `f64`.
let thresholds: [(f64, f64); 7] = [
(72.5, 1.0e-3),
(123.435, 1.5e-3),
(160.0, 7.1e-5),
@@ -169,7 +169,7 @@ fn test_the_model_jumps_at_the_species_threshold_altitudes() {
}
// Above 500 km there are none at all, which is the range most of a mission sits in.
let mut z = 500.0;
let mut z: f64 = 500.0;
while z < 1000.0 {
let step = 1.0e-9;
let below = model.gtd7d(&input(z - step)).total_mass;