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}