xad-rs 0.8.1

Exact automatic differentiation for Rust — forward-mode, reverse-mode, first- and second-order, with named variable support and a unified `Real` trait for mode-agnostic numerical code
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131

//! The `RealStats` extension trait — statistical functions (`erf`,
//! `erfc`, `norm_cdf`, `inv_norm_cdf`) for the active-scalar surface.
//!
//! Layered on top of [`crate::Real`] so callers that don't need
//! statistics can keep their bound at `R: Real`. Callers that price
//! options-style products write `R: Real + RealStats` (or just
//! `R: RealStats`, since `RealStats: Real`).
//!
//! # v0.6.0 scope
//!
//! Implemented for: `f64`, [`crate::AReal<f64>`], [`crate::Jet1<f64>`].
//!
//! Deliberately NOT implemented for `Jet2<f64>`, `Jet1Vec`, `Jet2Vec`
//! in v0.6.0 — those types lack inherent `erf`/`erfc`/`norm_cdf`/
//! `inv_norm_cdf` methods today, and deriving the second-order
//! chain-rule formulas (and the Newton-step derivative for
//! `inv_norm_cdf`) is substantial math work tracked under a separate
//! follow-up change `add-realstats-jet2-and-vec`.
//!
//! Numerical accuracy matches the underlying [`crate::math`] functions:
//! ~1.5e-7 for `erf` / `norm_cdf` (A&S 7.1.26); `inv_norm_cdf` uses
//! Acklam's rational approximation with comparable accuracy.

use crate::forward::Jet1;
use crate::real::Real;
use crate::reverse::AReal;

/// Statistical extension trait for active-scalar types. See module
/// docs for v0.6.0 scope.
pub trait RealStats: Real {
    /// Gaussian error function `erf(x) = (2/√π) ∫₀ˣ exp(-t²) dt`.
    fn erf(&self) -> Self;

    /// Complementary error function `erfc(x) = 1 − erf(x)`.
    fn erfc(&self) -> Self;

    /// Standard normal CDF `Φ(x) = ½ (1 + erf(x/√2))`.
    fn norm_cdf(&self) -> Self;

    /// Inverse standard normal CDF `Φ⁻¹(p)`.
    fn inv_norm_cdf(&self) -> Self;
}

// ----------------------------------------------------------------------------
// impl RealStats for f64
// ----------------------------------------------------------------------------

impl RealStats for f64 {
    #[inline]
    fn erf(&self) -> Self {
        crate::math::erf(*self)
    }
    #[inline]
    fn erfc(&self) -> Self {
        1.0 - crate::math::erf(*self)
    }
    #[inline]
    fn norm_cdf(&self) -> Self {
        crate::math::norm_cdf(*self)
    }
    #[inline]
    fn inv_norm_cdf(&self) -> Self {
        crate::math::inv_norm_cdf(*self)
    }
}

// ----------------------------------------------------------------------------
// impl RealStats for AReal<f64>
// ----------------------------------------------------------------------------

impl RealStats for AReal<f64> {
    #[inline]
    fn erf(&self) -> Self {
        crate::math::ad::erf(self)
    }
    #[inline]
    fn erfc(&self) -> Self {
        crate::math::ad::erfc(self)
    }
    #[inline]
    fn norm_cdf(&self) -> Self {
        crate::math::ad::norm_cdf(self)
    }
    #[inline]
    fn inv_norm_cdf(&self) -> Self {
        crate::math::ad::inv_norm_cdf(self)
    }
}

// ----------------------------------------------------------------------------
// impl RealStats for Jet1<f64>
// ----------------------------------------------------------------------------

impl RealStats for Jet1<f64> {
    #[inline]
    fn erf(&self) -> Self {
        crate::math::fwd::erf(self)
    }
    #[inline]
    fn erfc(&self) -> Self {
        crate::math::fwd::erfc(self)
    }
    #[inline]
    fn norm_cdf(&self) -> Self {
        crate::math::fwd::norm_cdf(self)
    }
    #[inline]
    fn inv_norm_cdf(&self) -> Self {
        crate::math::fwd::inv_norm_cdf(self)
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn f64_norm_cdf_at_zero_is_half() {
        assert!((<f64 as RealStats>::norm_cdf(&0.0) - 0.5).abs() < 1e-7);
    }

    #[test]
    fn f64_erf_at_zero_is_zero() {
        assert!(<f64 as RealStats>::erf(&0.0).abs() < 1e-7);
    }

    #[test]
    fn f64_erfc_at_zero_is_one() {
        assert!((<f64 as RealStats>::erfc(&0.0) - 1.0).abs() < 1e-7);
    }
}