num_dual/datatypes/
dual2_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 second order dual number for the calculation of Hessians.
12#[derive(Clone, Debug)]
13pub struct Dual2Vec<T: Scalar, D: Dim>
14where
15    DefaultAllocator: Allocator<U1, D> + Allocator<D, D>,
16{
17    /// Real part of the second order dual number
18    pub re: T,
19    /// Gradient part of the second order dual number
20    pub v1: Derivative<T, U1, D>,
21    /// Hessian part of the second order dual number
22    pub v2: Derivative<T, D, D>,
23}
24
25impl<T: Scalar + Copy, const N: usize> Copy for Dual2Vec<T, Const<N>> {}
26
27#[cfg(feature = "ndarray")]
28impl<T: DualNum, D: Dim> ndarray::ScalarOperand for Dual2Vec<T, D> where
29    DefaultAllocator: Allocator<U1, D> + Allocator<D, D>
30{
31}
32
33pub type Dual2SVec<T, const N: usize> = Dual2Vec<T, Const<N>>;
34pub type Dual2DVec<T> = Dual2Vec<T, Dyn>;
35pub type Dual2Vec32<D> = Dual2Vec<f32, D>;
36pub type Dual2Vec64<D> = Dual2Vec<f64, D>;
37pub type Dual2SVec32<const N: usize> = Dual2Vec<f32, Const<N>>;
38pub type Dual2SVec64<const N: usize> = Dual2Vec<f64, Const<N>>;
39pub type Dual2DVec32 = Dual2Vec<f32, Dyn>;
40pub type Dual2DVec64 = Dual2Vec<f64, Dyn>;
41
42impl<T: DualNum, D: Dim> Dual2Vec<T, D>
43where
44    DefaultAllocator: Allocator<U1, D> + Allocator<D, D>,
45{
46    /// Create a new second order dual number from its fields.
47    #[inline]
48    pub fn new(re: T, v1: Derivative<T, U1, D>, v2: Derivative<T, D, D>) -> Self {
49        Self { re, v1, v2 }
50    }
51}
52
53impl<T: DualNum, const N: usize> Dual2SVec<T, N> {
54    /// Set the derivative part of variable `index` to 1.
55    ///
56    /// For most cases, the [`hessian`](crate::hessian) function provides a convenient
57    /// interface to calculate derivatives. This function exists for the more edge cases
58    /// where more control over the variables is required.
59    /// ```
60    /// # use num_dual::Dual2SVec64;
61    /// # use nalgebra::{U1, U2, matrix};
62    /// let x: Dual2SVec64<2> = Dual2SVec64::from_re(5.0).derivative(0);
63    /// let y: Dual2SVec64<2> = Dual2SVec64::from_re(3.0).derivative(1);
64    /// let z = x * x * y;
65    /// assert_eq!(z.re, 75.0);                                                 // x²y
66    /// assert_eq!(z.v1.unwrap_generic(U1, U2), matrix![30.0, 25.0]);           // [2xy, x²]
67    /// assert_eq!(z.v2.unwrap_generic(U2, U2), matrix![6.0, 10.0; 10.0, 0.0]); // [2y, 2x; 2x, 0]
68    /// ```
69    #[inline]
70    pub fn derivative(mut self, index: usize) -> Self {
71        self.v1 = Derivative::derivative_generic(U1, Const::<N>, index);
72        self
73    }
74}
75
76impl<T: DualNum> Dual2DVec<T> {
77    /// Set the derivative part of variable `index` to 1.
78    ///
79    /// For most cases, the [`hessian`](crate::hessian) function provides a convenient interface
80    /// to calculate derivatives. This function exists for the more edge cases
81    /// where more control over the variables is required.
82    /// ```
83    /// # use num_dual::Dual2DVec64;
84    /// # use nalgebra::{Dyn, U1, dmatrix};
85    /// let x: Dual2DVec64 = Dual2DVec64::from_re(5.0).derivative(2, 0);
86    /// let y: Dual2DVec64 = Dual2DVec64::from_re(3.0).derivative(2, 1);
87    /// let z = &x * &x * y;
88    /// assert_eq!(z.re, 75.0);                                                          // x²y
89    /// assert_eq!(z.v1.unwrap_generic(U1, Dyn(2)), dmatrix![30.0, 25.0]);               // [2xy, x²]
90    /// assert_eq!(z.v2.unwrap_generic(Dyn(2), Dyn(2)), dmatrix![6.0, 10.0; 10.0, 0.0]); // [2y, 2x; 2x, 0]
91    /// ```
92    #[inline]
93    pub fn derivative(mut self, variables: usize, index: usize) -> Self {
94        self.v1 = Derivative::derivative_generic(U1, Dyn(variables), index);
95        self
96    }
97}
98
99impl<T: DualNum, D: Dim> Dual2Vec<T, D>
100where
101    DefaultAllocator: Allocator<U1, D> + Allocator<D, D>,
102{
103    /// Create a new second order dual number from the real part.
104    #[inline]
105    pub fn from_re(re: T) -> Self {
106        Self::new(re, Derivative::none(), Derivative::none())
107    }
108}
109
110/* chain rule */
111impl<T: DualNum, D: Dim> Dual2Vec<T, D>
112where
113    DefaultAllocator: Allocator<U1, D> + Allocator<D, D>,
114{
115    #[inline]
116    fn chain_rule(&self, f0: T, f1: T, f2: T) -> Self {
117        Self::new(
118            f0,
119            &self.v1 * f1.clone(),
120            &self.v2 * f1 + self.v1.tr_mul(&self.v1) * f2,
121        )
122    }
123}
124
125/* product rule */
126impl<T: DualNum, D: Dim> Mul<&Dual2Vec<T, D>> for &Dual2Vec<T, D>
127where
128    DefaultAllocator: Allocator<U1, D> + Allocator<D, D>,
129{
130    type Output = Dual2Vec<T, D>;
131    #[inline]
132    fn mul(self, other: &Dual2Vec<T, D>) -> Dual2Vec<T, D> {
133        Dual2Vec::new(
134            self.re.clone() * other.re.clone(),
135            &other.v1 * self.re.clone() + &self.v1 * other.re.clone(),
136            &other.v2 * self.re.clone()
137                + self.v1.tr_mul(&other.v1)
138                + other.v1.tr_mul(&self.v1)
139                + &self.v2 * other.re.clone(),
140        )
141    }
142}
143
144/* quotient rule */
145impl<T: DualNum, D: Dim> Div<&Dual2Vec<T, D>> for &Dual2Vec<T, D>
146where
147    DefaultAllocator: Allocator<U1, D> + Allocator<D, D>,
148{
149    type Output = Dual2Vec<T, D>;
150    #[inline]
151    fn div(self, other: &Dual2Vec<T, D>) -> Dual2Vec<T, D> {
152        let inv = other.re.recip();
153        let inv2 = inv.clone() * inv.clone();
154        Dual2Vec::new(
155            self.re.clone() * inv.clone(),
156            (&self.v1 * other.re.clone() - &other.v1 * self.re.clone()) * inv2.clone(),
157            &self.v2 * inv.clone()
158                - (&other.v2 * self.re.clone()
159                    + self.v1.tr_mul(&other.v1)
160                    + other.v1.tr_mul(&self.v1))
161                    * inv2.clone()
162                + other.v1.tr_mul(&other.v1)
163                    * ((T::one() + T::one()) * self.re.clone() * inv2 * inv),
164        )
165    }
166}
167
168/* string conversions */
169impl<T: DualNum, D: Dim> fmt::Display for Dual2Vec<T, D>
170where
171    DefaultAllocator: Allocator<U1, D> + Allocator<D, D>,
172{
173    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
174        write!(f, "{}", self.re)?;
175        self.v1.fmt(f, "ε1")?;
176        self.v2.fmt(f, "ε1²")
177    }
178}
179
180impl_second_derivatives!(Dual2Vec, [v1, v2], [D], [U1, D], [D, D]);
181impl_dual!(Dual2Vec, [v1, v2], [D], [U1, D], [D, D]);
182impl_nalgebra!(Dual2Vec, [v1, v2], [D], [U1, D], [D, D]);
num_dual/datatypes/dual2_vec.rs

num_dual/datatypes/
dual2_vec.rs
