num_dual/datatypes/
dual_vec.rs

1use crate::{Derivative, DualNum, DualNumFloat, DualStruct};
2use nalgebra::allocator::Allocator;
3use nalgebra::*;
4use num_traits::{FloatConst, FromPrimitive, Inv, Num, One, Signed, Zero};
5use std::fmt;
6use std::iter::{Product, Sum};
7use std::ops::{
8    Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Rem, RemAssign, Sub, SubAssign,
9};
10
11/// A vector dual number for the calculations of gradients or Jacobians.
12#[derive(Clone, Debug)]
13pub struct DualVec<T: Scalar, D: Dim>
14where
15    DefaultAllocator: Allocator<D>,
16{
17    /// Real part of the dual number
18    pub re: T,
19    /// Derivative part of the dual number
20    pub eps: Derivative<T, D, U1>,
21}
22
23#[cfg(feature = "ndarray")]
24impl<T: Scalar, D: Dim> ndarray::ScalarOperand for DualVec<T, D> where DefaultAllocator: Allocator<D>
25{}
26
27impl<T: Scalar + Copy, const N: usize> Copy for DualVec<T, Const<N>> {}
28
29pub type DualSVec<D, const N: usize> = DualVec<D, Const<N>>;
30pub type DualDVec<D> = DualVec<D, Dyn>;
31pub type DualVec32<D> = DualVec<f32, D>;
32pub type DualVec64<D> = DualVec<f64, D>;
33pub type DualSVec32<const N: usize> = DualVec<f32, Const<N>>;
34pub type DualSVec64<const N: usize> = DualVec<f64, Const<N>>;
35pub type DualDVec32 = DualVec<f32, Dyn>;
36pub type DualDVec64 = DualVec<f64, Dyn>;
37
38impl<T: DualNum, D: Dim> DualVec<T, D>
39where
40    DefaultAllocator: Allocator<D>,
41{
42    /// Create a new dual number from its fields.
43    #[inline]
44    pub fn new(re: T, eps: Derivative<T, D, U1>) -> Self {
45        Self { re, eps }
46    }
47}
48
49impl<T: DualNum, const N: usize> DualSVec<T, N> {
50    /// Set the derivative part of variable `index` to 1.
51    ///
52    /// For most cases, the [`gradient`](crate::gradient) function provides a convenient interface
53    /// to calculate derivatives. This function exists for the more edge cases
54    /// where more control over the variables is required.
55    /// ```
56    /// # use num_dual::DualSVec64;
57    /// # use nalgebra::{U1, U2, vector};
58    /// let x: DualSVec64<2> = DualSVec64::from_re(5.0).derivative(0);
59    /// let y: DualSVec64<2> = DualSVec64::from_re(3.0).derivative(1);
60    /// let z = x * x * y;
61    /// assert_eq!(z.re, 75.0);                                           // x²y
62    /// assert_eq!(z.eps.unwrap_generic(U2, U1), vector![30.0, 25.0]);    // [2xy, x²]
63    /// ```
64    #[inline]
65    pub fn derivative(mut self, index: usize) -> Self {
66        self.eps = Derivative::derivative_generic(Const::<N>, U1, index);
67        self
68    }
69}
70
71impl<T: DualNum> DualDVec<T> {
72    /// Set the derivative part of variable `index` to 1.
73    ///
74    /// For most cases, the [`gradient`](crate::gradient) function provides a convenient interface
75    /// to calculate derivatives. This function exists for the more edge cases
76    /// where more control over the variables is required.
77    /// ```
78    /// # use num_dual::DualDVec64;
79    /// # use nalgebra::{Dyn, U1, dvector};
80    /// let x: DualDVec64 = DualDVec64::from_re(5.0).derivative(2, 0);
81    /// let y: DualDVec64 = DualDVec64::from_re(3.0).derivative(2, 1);
82    /// let z = &x * &x * y;
83    /// assert_eq!(z.re, 75.0);                                               // x²y
84    /// assert_eq!(z.eps.unwrap_generic(Dyn(2), U1), dvector![30.0, 25.0]);   // [2xy, x²]
85    /// ```
86    #[inline]
87    pub fn derivative(mut self, variables: usize, index: usize) -> Self {
88        self.eps = Derivative::derivative_generic(Dyn(variables), U1, index);
89        self
90    }
91}
92
93impl<T: DualNum, D: Dim> DualVec<T, D>
94where
95    DefaultAllocator: Allocator<D>,
96{
97    /// Create a new dual number from the real part.
98    #[inline]
99    pub fn from_re(re: T) -> Self {
100        Self::new(re, Derivative::none())
101    }
102}
103
104/* chain rule */
105impl<T: DualNum, D: Dim> DualVec<T, D>
106where
107    DefaultAllocator: Allocator<D>,
108{
109    #[inline]
110    fn chain_rule(&self, f0: T, f1: T) -> Self {
111        Self::new(f0, &self.eps * f1)
112    }
113}
114
115/* product rule */
116impl<T: DualNum, D: Dim> Mul<&DualVec<T, D>> for &DualVec<T, D>
117where
118    DefaultAllocator: Allocator<D>,
119{
120    type Output = DualVec<T, D>;
121    #[inline]
122    fn mul(self, other: &DualVec<T, D>) -> Self::Output {
123        DualVec::new(
124            self.re.clone() * other.re.clone(),
125            &self.eps * other.re.clone() + &other.eps * self.re.clone(),
126        )
127    }
128}
129
130/* quotient rule */
131impl<T: DualNum, D: Dim> Div<&DualVec<T, D>> for &DualVec<T, D>
132where
133    DefaultAllocator: Allocator<D>,
134{
135    type Output = DualVec<T, D>;
136    #[inline]
137    fn div(self, other: &DualVec<T, D>) -> DualVec<T, D> {
138        let inv = other.re.recip();
139        DualVec::new(
140            self.re.clone() * inv.clone(),
141            (&self.eps * other.re.clone() - &other.eps * self.re.clone()) * inv.clone() * inv,
142        )
143    }
144}
145
146/* string conversions */
147impl<T: DualNum, D: Dim> fmt::Display for DualVec<T, D>
148where
149    DefaultAllocator: Allocator<D>,
150{
151    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
152        write!(f, "{}", self.re)?;
153        self.eps.fmt(f, "ε")
154    }
155}
156
157impl_first_derivatives!(DualVec, [eps], [D], [D]);
158impl_dual!(DualVec, [eps], [D], [D]);
159impl_nalgebra!(DualVec, [eps], [D], [D]);
num_dual/datatypes/dual_vec.rs

num_dual/datatypes/
dual_vec.rs
