num_dual/datatypes/
hyperdual_vec.rs

1use crate::{Derivative, DualNum, DualNumFloat, DualStruct};
2use nalgebra::allocator::Allocator;
3use nalgebra::{Const, DefaultAllocator, Dim, Dyn, Scalar, U1};
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 hyper-dual number for the calculation of partial Hessians.
12#[derive(Clone, Debug)]
13pub struct HyperDualVec<T: Scalar, M: Dim, N: Dim>
14where
15    DefaultAllocator: Allocator<M> + Allocator<M, N> + Allocator<U1, N>,
16{
17    /// Real part of the hyper-dual number
18    pub re: T,
19    /// Gradient part of the hyper-dual number
20    pub eps1: Derivative<T, M, U1>,
21    /// Gradient part of the hyper-dual number
22    pub eps2: Derivative<T, U1, N>,
23    /// Partial Hessian part of the hyper-dual number
24    pub eps1eps2: Derivative<T, M, N>,
25}
26
27impl<T: Scalar + Copy, const M: usize, const N: usize> Copy
28    for HyperDualVec<T, Const<M>, Const<N>>
29{
30}
31
32#[cfg(feature = "ndarray")]
33impl<T: Scalar, M: Dim, N: Dim> ndarray::ScalarOperand for HyperDualVec<T, M, N> where
34    DefaultAllocator: Allocator<M> + Allocator<M, N> + Allocator<U1, N>
35{
36}
37
38pub type HyperDualSVec<T, const M: usize, const N: usize> = HyperDualVec<T, Const<M>, Const<N>>;
39pub type HyperDualDVec<T> = HyperDualVec<T, Dyn, Dyn>;
40pub type HyperDualVec32<M, N> = HyperDualVec<f32, M, N>;
41pub type HyperDualVec64<M, N> = HyperDualVec<f64, M, N>;
42pub type HyperDualSVec32<const M: usize, const N: usize> = HyperDualVec<f32, Const<M>, Const<N>>;
43pub type HyperDualSVec64<const M: usize, const N: usize> = HyperDualVec<f64, Const<M>, Const<N>>;
44pub type HyperDualDVec32 = HyperDualVec<f32, Dyn, Dyn>;
45pub type HyperDualDVec64 = HyperDualVec<f64, Dyn, Dyn>;
46
47impl<T: DualNum, M: Dim, N: Dim> HyperDualVec<T, M, N>
48where
49    DefaultAllocator: Allocator<M> + Allocator<M, N> + Allocator<U1, N>,
50{
51    /// Create a new hyper-dual number from its fields.
52    #[inline]
53    pub fn new(
54        re: T,
55        eps1: Derivative<T, M, U1>,
56        eps2: Derivative<T, U1, N>,
57        eps1eps2: Derivative<T, M, N>,
58    ) -> Self {
59        Self {
60            re,
61            eps1,
62            eps2,
63            eps1eps2,
64        }
65    }
66}
67
68impl<T: DualNum, M: Dim, N: Dim> HyperDualVec<T, M, N>
69where
70    DefaultAllocator: Allocator<M> + Allocator<M, N> + Allocator<U1, N>,
71{
72    /// Create a new hyper-dual number from the real part.
73    #[inline]
74    pub fn from_re(re: T) -> Self {
75        Self::new(
76            re,
77            Derivative::none(),
78            Derivative::none(),
79            Derivative::none(),
80        )
81    }
82}
83
84impl<T: DualNum, const M: usize, const N: usize> HyperDualSVec<T, M, N> {
85    /// Set the 1st dimension derivative of variable `index` to 1.
86    ///
87    /// For most cases, the [`partial_hessian`](crate::partial_hessian) function provides a
88    /// convenient interface to calculate derivatives. This function exists for the more edge
89    /// cases where more control over the variables is required.
90    #[inline]
91    pub fn derivative1(mut self, index: usize) -> Self {
92        self.eps1 = Derivative::derivative_generic(Const::<M>, U1, index);
93        self
94    }
95
96    /// Set the 2nd dimension derivative of variable `index` to 1.
97    ///
98    /// For most cases, the [`partial_hessian`](crate::partial_hessian) function provides a
99    /// convenient interface to calculate derivatives. This function exists for the more edge
100    /// cases where more control over the variables is required.
101    #[inline]
102    pub fn derivative2(mut self, index: usize) -> Self {
103        self.eps2 = Derivative::derivative_generic(U1, Const::<N>, index);
104        self
105    }
106}
107
108impl<T: DualNum> HyperDualDVec<T> {
109    /// Set the 1st dimension derivative part of variable `index` to 1.
110    ///
111    /// For most cases, the [`partial_hessian`](crate::partial_hessian) function provides a
112    /// convenient interface to calculate derivatives. This function exists for the more edge
113    /// cases where more control over the variables is required.
114    #[inline]
115    pub fn derivative1(mut self, variables: usize, index: usize) -> Self {
116        self.eps1 = Derivative::derivative_generic(Dyn(variables), U1, index);
117        self
118    }
119
120    /// Set the 2nd dimension derivative part of variable `index` to 1.
121    ///
122    /// For most cases, the [`partial_hessian`](crate::partial_hessian) function provides a
123    /// convenient interface to calculate derivatives. This function exists for the more edge
124    /// cases where more control over the variables is required.
125    #[inline]
126    pub fn derivative2(mut self, variables: usize, index: usize) -> Self {
127        self.eps2 = Derivative::derivative_generic(U1, Dyn(variables), index);
128        self
129    }
130}
131
132/* chain rule */
133impl<T: DualNum, M: Dim, N: Dim> HyperDualVec<T, M, N>
134where
135    DefaultAllocator: Allocator<M> + Allocator<M, N> + Allocator<U1, N>,
136{
137    #[inline]
138    fn chain_rule(&self, f0: T, f1: T, f2: T) -> Self {
139        Self::new(
140            f0,
141            &self.eps1 * f1.clone(),
142            &self.eps2 * f1.clone(),
143            &self.eps1eps2 * f1 + &self.eps1 * &self.eps2 * f2,
144        )
145    }
146}
147
148/* product rule */
149impl<T: DualNum, M: Dim, N: Dim> Mul<&HyperDualVec<T, M, N>> for &HyperDualVec<T, M, N>
150where
151    DefaultAllocator: Allocator<M> + Allocator<M, N> + Allocator<U1, N>,
152{
153    type Output = HyperDualVec<T, M, N>;
154    #[inline]
155    fn mul(self, other: &HyperDualVec<T, M, N>) -> HyperDualVec<T, M, N> {
156        HyperDualVec::new(
157            self.re.clone() * other.re.clone(),
158            &other.eps1 * self.re.clone() + &self.eps1 * other.re.clone(),
159            &other.eps2 * self.re.clone() + &self.eps2 * other.re.clone(),
160            &other.eps1eps2 * self.re.clone()
161                + &self.eps1 * &other.eps2
162                + &other.eps1 * &self.eps2
163                + &self.eps1eps2 * other.re.clone(),
164        )
165    }
166}
167
168/* quotient rule */
169impl<T: DualNum, M: Dim, N: Dim> Div<&HyperDualVec<T, M, N>> for &HyperDualVec<T, M, N>
170where
171    DefaultAllocator: Allocator<M> + Allocator<M, N> + Allocator<U1, N>,
172{
173    type Output = HyperDualVec<T, M, N>;
174    #[inline]
175    fn div(self, other: &HyperDualVec<T, M, N>) -> HyperDualVec<T, M, N> {
176        let inv = other.re.recip();
177        let inv2 = inv.clone() * &inv;
178        HyperDualVec::new(
179            self.re.clone() * &inv,
180            (&self.eps1 * other.re.clone() - &other.eps1 * self.re.clone()) * inv2.clone(),
181            (&self.eps2 * other.re.clone() - &other.eps2 * self.re.clone()) * inv2.clone(),
182            &self.eps1eps2 * inv.clone()
183                - (&other.eps1eps2 * self.re.clone()
184                    + &self.eps1 * &other.eps2
185                    + &other.eps1 * &self.eps2)
186                    * inv2.clone()
187                + &other.eps1
188                    * &other.eps2
189                    * ((T::one() + T::one()) * self.re.clone() * inv2 * inv),
190        )
191    }
192}
193
194/* string conversions */
195impl<T: DualNum, M: Dim, N: Dim> fmt::Display for HyperDualVec<T, M, N>
196where
197    DefaultAllocator: Allocator<M> + Allocator<M, N> + Allocator<U1, N>,
198{
199    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
200        write!(f, "{}", self.re)?;
201        self.eps1.fmt(f, "ε1")?;
202        self.eps2.fmt(f, "ε2")?;
203        self.eps1eps2.fmt(f, "ε1ε2")
204    }
205}
206
207impl_second_derivatives!(
208    HyperDualVec,
209    [eps1, eps2, eps1eps2],
210    [M, N],
211    [M],
212    [M, N],
213    [U1, N]
214);
215impl_dual!(
216    HyperDualVec,
217    [eps1, eps2, eps1eps2],
218    [M, N],
219    [M],
220    [M, N],
221    [U1, N]
222);
223impl_nalgebra!(
224    HyperDualVec,
225    [eps1, eps2, eps1eps2],
226    [M, N],
227    [M],
228    [M, N],
229    [U1, N]
230);
num_dual/datatypes/hyperdual_vec.rs

num_dual/datatypes/
hyperdual_vec.rs
