raddy/
lib.rs

1#![doc = include_str!("../README.md")]
2extern crate nalgebra as na;
3
4/// Comparison operations and utilities for AD values.
5pub mod compare;
6
7/// Factory functions for creating AD values and vectors.
8pub mod make;
9
10/// Miscellaneous utilities and experimental features.
11mod misc;
12
13/// Scalar operations, operator traits, and field implementations.
14/// Please Note that all `unimplemented!` methods are not intended for use.
15/// If any operation encountered these, please raise an issue.
16pub mod scalar;
17
18/// Sparse matrix differentiation functionalities.
19pub mod sparse;
20
21/// Unit tests and validation code.
22#[cfg(test)]
23mod test;
24
25/// Core type definitions and AD value implementations.
26pub mod types;
27
28// ######################################################################################
29// ################################### Implementation ###################################
30// ######################################################################################
31
32use na::{SMatrix, SVector};
33use types::{mat, vec};
34
35// ################################### Data Structure ###################################
36
37/// Automatic differentiation value tracking first and second derivatives
38///
39/// # Value getters:
40/// - `value() -> f64`: Returns the current numerical value
41/// - `grad() -> SVector<f64, N>`: Returns the gradient vector
42/// - `hess() -> SMatrix<f64, N, N>`: Returns the Hessian matrix
43///
44/// # Type Parameters
45/// * `N` - The dimension of the input space (number of variables)
46///
47/// # Fields (private)
48/// * `value` - The current value of the function
49/// * `grad` - The gradient (first derivatives) as a vector
50/// * `hess` - The Hessian matrix (second derivatives)
51#[derive(Debug, Clone)]
52pub struct Ad<const N: usize> {
53    pub(crate) value: f64,
54    pub(crate) grad: vec<N>,
55    pub(crate) hess: mat<N>,
56}
57
58// ################################### Accessors ###################################
59
60impl<const N: usize> Ad<N> {
61    /// Returns the current value of the AD variable
62    pub fn value(&self) -> f64 {
63        self.value
64    }
65
66    /// Returns the gradient (first derivatives) of the AD variable
67    ///
68    /// # Returns
69    /// The gradient (A vector containing the partial derivatives with respect to each input variable).
70    pub fn grad(&self) -> vec<N> {
71        self.grad.clone()
72    }
73
74    /// Returns the Hessian matrix (second derivatives) of the AD variable
75    ///
76    /// # Returns
77    /// The [hessian](https://en.wikipedia.org/wiki/Hessian_matrix).
78    pub fn hess(&self) -> mat<N> {
79        self.hess.clone()
80    }
81}
82
83/// Trait for extracting numerical values from AD matrices
84///
85/// # Type Parameters
86/// * `R` - Number of rows
87/// * `C` - Number of columns
88pub trait GetValue<const R: usize, const C: usize> {
89    /// The type of the extracted numerical values
90    type Value;
91
92    /// Extracts the numerical values from an AD matrix
93    fn value(&self) -> Self::Value;
94}
95
96impl<const N: usize, const R: usize, const C: usize> GetValue<R, C> for SMatrix<Ad<N>, R, C> {
97    type Value = SMatrix<f64, R, C>;
98    fn value(&self) -> Self::Value {
99        let mut val = Self::Value::zeros();
100        for r in 0..R {
101            for c in 0..C {
102                val[(r, c)] = self[(r, c)].value;
103            }
104        }
105        val
106    }
107}
108
109// ################################### Public Constructors ###################################
110
111impl Ad<1> {
112    /// Creates an AD scalar with explicitly specified value, gradient and Hessian
113    ///
114    /// # Arguments
115    /// * `value` - The scalar value
116    /// * `grad` - The gradient (first derivative)
117    /// * `hess` - The Hessian (second derivative)
118    ///
119    /// # Returns
120    /// A new Ad<1> instance with the specified properties
121    pub fn given_scalar(value: f64, grad: f64, hess: f64) -> Self {
122        Self {
123            value,
124            grad: vec::from_row_slice(&[grad]),
125            hess: mat::from_row_slice(&[hess]),
126        }
127    }
128
129    /// Creates an active scalar AD value with unit gradient
130    ///
131    /// # Arguments
132    /// * `value` - The scalar value
133    ///
134    /// # Returns
135    /// A new Ad<1> instance that is active (gradient = 1.0)
136    pub fn active_scalar(value: f64) -> Self {
137        let mut res = Self::_zeroed();
138
139        res.value = value;
140        res.grad[0] = 1.0;
141
142        res
143    }
144}
145
146impl<const N: usize> Ad<N> {
147    /// Creates an inactive scalar AD value with zero gradient and Hessian
148    ///
149    /// # Arguments
150    /// * `value` - The scalar value
151    ///
152    /// # Returns
153    /// A new `Ad<N>` instance that is inactive (gradient = 0)
154    pub fn inactive_scalar(value: f64) -> Self {
155        let mut res = Self::_zeroed();
156        res.value = value;
157        res
158    }
159
160    /// Creates a vector of inactive AD values from a vector of f64 values
161    ///
162    /// # Arguments
163    /// * `values` - Input vector of numerical values
164    ///
165    /// # Type Parameters
166    /// * `L` - Length of the output vector
167    ///
168    /// # Returns
169    /// A vector of inactive AD values
170    pub fn inactive_vector<const L: usize>(values: &SVector<f64, N>) -> SVector<Self, L> {
171        let mut scalars = Vec::new();
172        for i in 0..N {
173            let scalar = Self::inactive_scalar(values[i]);
174            scalars.push(scalar);
175        }
176
177        SVector::from_column_slice(&scalars)
178    }
179
180    /// Creates a vector of inactive AD values from a slice of f64 values
181    ///
182    /// # Arguments
183    /// * `values` - Slice of numerical values
184    ///
185    /// # Type Parameters
186    /// * `L` - Length of the output vector
187    ///
188    /// # Returns
189    /// A vector of inactive AD values
190    ///
191    /// # Panics
192    /// If the slice length doesn't match the input dimension N
193    pub fn inactive_from_slice<const L: usize>(values: &[f64]) -> SVector<Self, L> {
194        assert_eq!(
195            values.len(),
196            N,
197            "Slice length mismatch: expected {}, got {}",
198            N,
199            values.len()
200        );
201        Self::inactive_vector(&SVector::from_column_slice(values))
202    }
203
204    /// Creates an AD value with explicitly specified value, gradient and Hessian
205    ///
206    /// # Arguments
207    /// * `value` - The scalar value
208    /// * `grad` - The gradient vector
209    /// * `hess` - The Hessian matrix
210    ///
211    /// # Returns
212    /// A new `Ad<N>` instance with the specified properties
213    pub fn given_vector(value: f64, grad: &vec<N>, hess: &mat<N>) -> Self {
214        Self {
215            value,
216            grad: grad.clone(),
217            hess: hess.clone(),
218        }
219    }
220
221    /// Creates a vector of active AD values from a vector of f64 values
222    ///
223    /// # Arguments
224    /// * `values` - Input vector of numerical values
225    ///
226    /// # Returns
227    /// A vector of active AD values where each element has unit gradient
228    /// in its corresponding dimension
229    pub fn active_vector(vector: &SVector<f64, N>) -> SVector<Self, N> {
230        let mut scalars = Vec::new();
231        for i in 0..N {
232            let scalar = Self::_active_scalar_with_index(vector[i], i);
233            scalars.push(scalar);
234        }
235
236        SVector::from_column_slice(&scalars)
237    }
238
239    /// Creates a vector of active AD values from a slice of f64 values
240    ///
241    /// # Arguments
242    /// * `values` - Slice of numerical values
243    ///
244    /// # Returns
245    /// A vector of active AD values
246    ///
247    /// # Panics
248    /// If the slice length doesn't match the input dimension N
249    pub fn active_from_slice(values: &[f64]) -> SVector<Self, N> {
250        assert_eq!(
251            values.len(),
252            N,
253            "Slice length mismatch: expected {}, got {}",
254            N,
255            values.len()
256        );
257        Self::active_vector(&SVector::from_column_slice(values))
258    }
259}
260
261// ################################### Private Constructors ###################################
262
263impl<const N: usize> Ad<N> {
264    fn _active_scalar_with_index(value: f64, index: usize) -> Self {
265        let mut res = Self::_zeroed();
266
267        res.value = value;
268        res.grad[index] = 1.0;
269
270        res
271    }
272
273    fn _zeroed() -> Self {
274        Self {
275            value: 0.0,
276            grad: vec::zeros(),
277            hess: mat::zeros(),
278        }
279    }
280}
281
282// ################################### Utils ###################################
283
284impl<const N: usize> Ad<N> {
285    fn chain(
286        value: f64, // f
287        d: f64,     // df/da
288        d2: f64,    // ddf/daa
289        a: &Self,
290    ) -> Self {
291        let mut res = Self::_zeroed();
292
293        res.value = value;
294        res.grad = d * &a.grad;
295        res.hess = d2 * &a.grad * &a.grad.transpose() + d * a.hess;
296
297        res
298    }
299}