From b3d95874b03bffb31f5edb75978487f6f34f2554 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Thu, 5 Sep 2024 20:11:36 +0200 Subject: [PATCH 1/2] replace allow with expect --- src/derivative.rs | 2 +- src/dual2_vec.rs | 2 +- src/dual_size.rs | 309 +++++++++++++++++++++++++++++++++++ src/dual_vec.rs | 2 +- src/hyperdual_vec.rs | 4 +- src/hyperhyperdual.rs | 6 +- src/lib.rs | 2 +- src/macros.rs | 16 -- src/python/hyperdual.rs | 1 - src/python/hyperhyperdual.rs | 6 +- 10 files changed, 321 insertions(+), 29 deletions(-) create mode 100644 src/dual_size.rs diff --git a/src/derivative.rs b/src/derivative.rs index c878ebd..d9513dd 100644 --- a/src/derivative.rs +++ b/src/derivative.rs @@ -140,7 +140,7 @@ where } impl, F> Derivative { - #[allow(clippy::self_named_constructors)] + #[expect(clippy::self_named_constructors)] pub fn derivative() -> Self { Self::some(SVector::identity()) } diff --git a/src/dual2_vec.rs b/src/dual2_vec.rs index 7fc3db7..8428f32 100644 --- a/src/dual2_vec.rs +++ b/src/dual2_vec.rs @@ -89,7 +89,7 @@ where } /// Variant of [hessian] for fallible functions. -#[allow(clippy::type_complexity)] +#[expect(clippy::type_complexity)] pub fn try_hessian, F: DualNumFloat, E, D: Dim>( g: G, x: OVector, diff --git a/src/dual_size.rs b/src/dual_size.rs new file mode 100644 index 0000000..2fe0bee --- /dev/null +++ b/src/dual_size.rs @@ -0,0 +1,309 @@ +use crate::{DualNum, DualNumFloat, OneZero, Scale}; +use nalgebra::allocator::Allocator; +use nalgebra::*; +use num_traits::{Float, Inv, One}; +use std::convert::Infallible; +use std::fmt; +use std::iter::{Product, Sum}; +use std::marker::PhantomData; +use std::ops::{ + Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Rem, RemAssign, Sub, SubAssign, +}; + +/// A dual number for the calculations of gradients or Jacobians. +#[derive(PartialEq, Eq, Clone, Debug)] +pub struct DualVec, F, N: Dim> +where + DefaultAllocator: Allocator, +{ + /// Real part of the dual number + pub re: T, + /// Derivative part of the dual number + pub eps: OVector, + f: PhantomData, +} + +impl + Copy, F: Copy, const N: usize> Copy for DualVec> {} + +pub type DualVec32 = DualVec; +pub type DualVec64 = DualVec; +pub type DualSVec32 = DualVec>; +pub type DualSVec64 = DualVec>; +pub type DualDVec32 = DualVec; +pub type DualDVec64 = DualVec; +pub type Dual = DualVec; +pub type Dual32 = Dual; +pub type Dual64 = Dual; + +impl, F, N: Dim> DualVec +where + DefaultAllocator: Allocator, +{ + /// Create a new dual number from its fields. + #[inline] + pub fn new(re: T, eps: OVector) -> Self { + Self { + re, + eps, + f: PhantomData, + } + } +} + +impl, F, N: Dim> OneZero for DualVec +where + DefaultAllocator: Allocator, +{ + fn zero(&self) -> Self { + let (n, _) = self.eps.shape_generic(); + Self::new( + self.re.zero(), + OVector::from_element_generic(n, U1, self.re.zero()), + ) + } + + fn one(&self) -> Self { + let (n, _) = self.eps.shape_generic(); + Self::new( + self.re.one(), + OVector::from_element_generic(n, U1, self.re.zero()), + ) + } +} + +impl, F> Dual { + /// Create a new scalar dual number from its fields. + #[inline] + pub fn new_scalar(re: T, eps: T) -> Self { + Self::new(re, SVector::from_element(eps)) + } +} + +// impl + Zero, F, N: Dim> DualVec +// where +// DefaultAllocator: Allocator, +// { +// /// Create a new dual number from the real part. +// #[inline] +// pub fn from_re(re: T) -> Self { +// Self::new(re, OVector::zeros_generic()) +// } +// } + +impl + One, F> Dual { + /// Set the derivative part to 1. + /// ``` + /// # use num_dual::{Dual64, DualNum}; + /// let x = Dual64::from_re(5.0).derivative().powi(2); + /// assert_eq!(x.re, 25.0); + /// assert_eq!(x.eps.unwrap(), 10.0); + /// ``` + #[inline] + pub fn derivative(x: T) -> Self { + let one = x.one(); + Self::new(x, SVector::from_element(one)) + } +} + +/// Calculate the first derivative of a scalar function. +/// ``` +/// # use num_dual::{first_derivative, DualNum}; +/// let (f, df) = first_derivative(|x| x.powi(2), 5.0); +/// assert_eq!(f, 25.0); +/// assert_eq!(df, 10.0); +/// ``` +pub fn first_derivative, F>(g: G, x: T) -> (T, T) +where + G: FnOnce(Dual) -> Dual, +{ + try_first_derivative(|x| Ok::<_, Infallible>(g(x)), x).unwrap() +} + +/// Variant of [first_derivative] for fallible functions. +pub fn try_first_derivative, F, E>(g: G, x: T) -> Result<(T, T), E> +where + G: FnOnce(Dual) -> Result, E>, +{ + let one = x.one(); + let x = Dual::new_scalar(x, one); + g(x).map(|r| (r.re, r.eps.data.0[0][0].clone())) +} + +/// Calculate the gradient of a scalar function +/// ``` +/// # use approx::assert_relative_eq; +/// # use num_dual::{gradient, DualNum, DualSVec64}; +/// # use nalgebra::SVector; +/// let v = SVector::from([4.0, 3.0]); +/// let fun = |v: SVector, 2>| (v[0].powi(2) + v[1].powi(2)).sqrt(); +/// let (f, g) = gradient(fun, v); +/// assert_eq!(f, 5.0); +/// assert_relative_eq!(g[0], 0.8); +/// assert_relative_eq!(g[1], 0.6); +/// ``` +/// +/// The variable vector can also be dynamically sized +/// ``` +/// # use approx::assert_relative_eq; +/// # use num_dual::{gradient, DualNum, DualDVec64}; +/// # use nalgebra::DVector; +/// let v = DVector::repeat(4, 2.0); +/// let fun = |v: DVector| v.iter().map(|v| v * v).sum::().sqrt(); +/// let (f, g) = gradient(fun, v); +/// assert_eq!(f, 4.0); +/// assert_relative_eq!(g[0], 0.5); +/// assert_relative_eq!(g[1], 0.5); +/// assert_relative_eq!(g[2], 0.5); +/// assert_relative_eq!(g[3], 0.5); +/// ``` +pub fn gradient, F: DualNumFloat, N: Dim>( + g: G, + x: OVector, +) -> (T, OVector) +where + G: FnOnce(OVector, N>) -> DualVec, + DefaultAllocator: Allocator + Allocator, N>, +{ + try_gradient(|x| Ok::<_, Infallible>(g(x)), x).unwrap() +} + +/// Variant of [gradient] for fallible functions. +pub fn try_gradient, F: DualNumFloat, E, N: Dim>( + g: G, + x: OVector, +) -> Result<(T, OVector), E> +where + G: FnOnce(OVector, N>) -> Result, E>, + DefaultAllocator: Allocator + Allocator, N>, +{ + let (n, _) = x.shape_generic(); + let x = x.map_with_location(|i, _, x| { + let mut eps = OVector::from_element_generic(n, U1, x.zero()); + eps[i] = x.one(); + DualVec::new(x, eps) + }); + g(x).map(|res| (res.re, res.eps)) +} + +/// Calculate the Jacobian of a vector function. +/// ``` +/// # use num_dual::{jacobian, DualSVec64, DualNum}; +/// # use nalgebra::SVector; +/// let xy = SVector::from([5.0, 3.0, 2.0]); +/// let fun = |xy: SVector, 3>| SVector::from([ +/// xy[0] * xy[1].powi(3) * xy[2], +/// xy[0].powi(2) * xy[1] * xy[2].powi(2) +/// ]); +/// let (f, jac) = jacobian(fun, xy); +/// assert_eq!(f[0], 270.0); // xy³z +/// assert_eq!(f[1], 300.0); // x²yz² +/// assert_eq!(jac[(0,0)], 54.0); // y³z +/// assert_eq!(jac[(0,1)], 270.0); // 3xy²z +/// assert_eq!(jac[(0,2)], 135.0); // xy³ +/// assert_eq!(jac[(1,0)], 120.0); // 2xyz² +/// assert_eq!(jac[(1,1)], 100.0); // x²z² +/// assert_eq!(jac[(1,2)], 300.0); // 2x²yz +/// ``` +pub fn jacobian, F: DualNumFloat, M: Dim, N: Dim>( + g: G, + x: OVector, +) -> (OVector, OMatrix) +where + G: FnOnce(OVector, N>) -> OVector, M>, + DefaultAllocator: Allocator, M> + + Allocator + + Allocator + + Allocator + + Allocator, N> + + Allocator, N> + + Allocator, M>, +{ + try_jacobian(|x| Ok::<_, Infallible>(g(x)), x).unwrap() +} + +/// Variant of [jacobian] for fallible functions. +#[allow(clippy::type_complexity)] +pub fn try_jacobian, F: DualNumFloat, E, M: Dim, N: Dim>( + g: G, + x: OVector, +) -> Result<(OVector, OMatrix), E> +where + G: FnOnce(OVector, N>) -> Result, M>, E>, + DefaultAllocator: Allocator, M> + + Allocator + + Allocator + + Allocator + + Allocator, N> + + Allocator, N> + + Allocator, M>, +{ + let (n, _) = x.shape_generic(); + let x = x.map_with_location(|i, _, x| { + let mut eps = OVector::from_element_generic(n, U1, x.zero()); + eps[i] = x.one(); + DualVec::new(x, eps) + }); + g(x).map(|res| { + let eps = OMatrix::from_rows(res.map(|res| res.eps.transpose()).as_slice()); + (res.map(|r| r.re), eps) + }) +} + +/* chain rule */ +impl, F: Float, N: Dim> DualVec +where + DefaultAllocator: Allocator, +{ + #[inline] + fn chain_rule(&self, f0: T, f1: T) -> Self { + Self::new(f0, self.eps.clone() * f1) + } +} + +/* product rule */ +impl<'a, 'b, T: DualNum, F: Float, N: Dim> Mul<&'a DualVec> for &'b DualVec +where + DefaultAllocator: Allocator, +{ + type Output = DualVec; + #[inline] + fn mul(self, other: &DualVec) -> Self::Output { + DualVec::new( + self.re.clone() * other.re.clone(), + &self.eps * other.re.clone() + &other.eps * self.re.clone(), + ) + } +} + +/* quotient rule */ +impl<'a, 'b, T: DualNum, F: Float, N: Dim> Div<&'a DualVec> for &'b DualVec +where + DefaultAllocator: Allocator, +{ + type Output = DualVec; + #[inline] + fn div(self, other: &DualVec) -> DualVec { + let inv = other.re.recip(); + DualVec::new( + self.re.clone() * inv.clone(), + (self.eps.clone() * other.re.clone() - other.eps.clone() * self.re.clone()) + * inv.clone() + * inv, + ) + } +} + +/* string conversions */ +impl, F, N: Dim> fmt::Display for DualVec +where + DefaultAllocator: Allocator, +{ + fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { + todo!() + // write!(f, "{}", self.re)?; + // self.eps.fmt(f, "ε") + } +} + +impl_first_derivatives2!(DualVec, [eps], [D]); +impl_dual2!(DualVec, [eps], [D]); diff --git a/src/dual_vec.rs b/src/dual_vec.rs index 265cf71..6c055b1 100644 --- a/src/dual_vec.rs +++ b/src/dual_vec.rs @@ -146,7 +146,7 @@ where } /// Variant of [jacobian] for fallible functions. -#[allow(clippy::type_complexity)] +#[expect(clippy::type_complexity)] pub fn try_jacobian, F: DualNumFloat, E, M: Dim, N: Dim>( g: G, x: OVector, diff --git a/src/hyperdual_vec.rs b/src/hyperdual_vec.rs index b2b38cd..45c3565 100644 --- a/src/hyperdual_vec.rs +++ b/src/hyperdual_vec.rs @@ -96,7 +96,7 @@ where /// assert_relative_eq!(d2fdxdy[0], -0.032); /// assert_relative_eq!(d2fdxdy[1], -0.024); /// ``` -#[allow(clippy::type_complexity)] +#[expect(clippy::type_complexity)] pub fn partial_hessian, F: DualNumFloat, M: Dim, N: Dim>( g: G, x: OVector, @@ -113,7 +113,7 @@ where } /// Variant of [partial_hessian] for fallible functions. -#[allow(clippy::type_complexity)] +#[expect(clippy::type_complexity)] pub fn try_partial_hessian, F: DualNumFloat, E, M: Dim, N: Dim>( g: G, x: OVector, diff --git a/src/hyperhyperdual.rs b/src/hyperhyperdual.rs index a0a3db5..0a64b9c 100644 --- a/src/hyperhyperdual.rs +++ b/src/hyperhyperdual.rs @@ -38,7 +38,7 @@ pub type HyperHyperDual64 = HyperHyperDual; impl, F> HyperHyperDual { /// Create a new hyper-hyper-dual number from its fields. #[inline] - #[allow(clippy::too_many_arguments)] + #[expect(clippy::too_many_arguments)] pub fn new( re: T, eps1: T, @@ -133,7 +133,7 @@ where } /// Variant of [third_partial_derivative] for fallible functions. -#[allow(clippy::type_complexity)] +#[expect(clippy::type_complexity)] pub fn try_third_partial_derivative, F, E>( g: G, x: T, @@ -199,7 +199,7 @@ where } /// Variant of [third_partial_derivative_vec] for fallible functions. -#[allow(clippy::type_complexity)] +#[expect(clippy::type_complexity)] pub fn try_third_partial_derivative_vec, F, E>( g: G, x: &[T], diff --git a/src/lib.rs b/src/lib.rs index bfd158e..cbd8108 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -40,7 +40,7 @@ //! ``` #![warn(clippy::all)] -#![allow(clippy::needless_range_loop)] +#![warn(clippy::allow_attributes)] use num_traits::{Float, FloatConst, FromPrimitive, Inv, NumAssignOps, NumOps, Signed}; use std::fmt; diff --git a/src/macros.rs b/src/macros.rs index aab4e46..ec1b6fe 100644 --- a/src/macros.rs +++ b/src/macros.rs @@ -557,82 +557,66 @@ macro_rules! impl_float_const { $($(DefaultAllocator: Allocator<$dim> + Allocator + Allocator<$dim, $dim>,)* DefaultAllocator: Allocator<$($dim,)*>)? { - #[allow(non_snake_case)] fn E() -> Self { Self::from(F::E()) } - #[allow(non_snake_case)] fn FRAC_1_PI() -> Self { Self::from(F::FRAC_1_PI()) } - #[allow(non_snake_case)] fn FRAC_1_SQRT_2() -> Self { Self::from(F::FRAC_1_SQRT_2()) } - #[allow(non_snake_case)] fn FRAC_2_PI() -> Self { Self::from(F::FRAC_2_PI()) } - #[allow(non_snake_case)] fn FRAC_2_SQRT_PI() -> Self { Self::from(F::FRAC_2_SQRT_PI()) } - #[allow(non_snake_case)] fn FRAC_PI_2() -> Self { Self::from(F::FRAC_PI_2()) } - #[allow(non_snake_case)] fn FRAC_PI_3() -> Self { Self::from(F::FRAC_PI_3()) } - #[allow(non_snake_case)] fn FRAC_PI_4() -> Self { Self::from(F::FRAC_PI_4()) } - #[allow(non_snake_case)] fn FRAC_PI_6() -> Self { Self::from(F::FRAC_PI_6()) } - #[allow(non_snake_case)] fn FRAC_PI_8() -> Self { Self::from(F::FRAC_PI_8()) } - #[allow(non_snake_case)] fn LN_10() -> Self { Self::from(F::LN_10()) } - #[allow(non_snake_case)] fn LN_2() -> Self { Self::from(F::LN_2()) } - #[allow(non_snake_case)] fn LOG10_E() -> Self { Self::from(F::LOG10_E()) } - #[allow(non_snake_case)] fn LOG2_E() -> Self { Self::from(F::LOG2_E()) } - #[allow(non_snake_case)] fn PI() -> Self { Self::from(F::PI()) } - #[allow(non_snake_case)] fn SQRT_2() -> Self { Self::from(F::SQRT_2()) } diff --git a/src/python/hyperdual.rs b/src/python/hyperdual.rs index 67c4eb8..70a148a 100644 --- a/src/python/hyperdual.rs +++ b/src/python/hyperdual.rs @@ -161,7 +161,6 @@ macro_rules! impl_partial_hessian { /// Returns /// ------- /// function value, gradient w.r.t. x, gradient w.r.t. y, and partial Hessian - #[allow(clippy::type_complexity)] pub fn partial_hessian( f: &Bound<'_, PyAny>, x: &Bound<'_, PyAny>, diff --git a/src/python/hyperhyperdual.rs b/src/python/hyperhyperdual.rs index 8a1f9af..8c0fe5e 100644 --- a/src/python/hyperhyperdual.rs +++ b/src/python/hyperhyperdual.rs @@ -9,7 +9,7 @@ use pyo3::prelude::*; pub struct PyHyperHyperDual64(HyperHyperDual64); #[pymethods] -#[allow(clippy::too_many_arguments)] +#[expect(clippy::too_many_arguments)] impl PyHyperHyperDual64 { #[new] fn new( @@ -77,7 +77,7 @@ impl_dual_num!(PyHyperHyperDual64, HyperHyperDual64, f64); /// second partial derivative w.r.t. x and z /// second partial derivative w.r.t. y and z /// third partial derivative -#[allow(clippy::type_complexity)] +#[expect(clippy::type_complexity)] pub fn third_partial_derivative( f: &Bound<'_, PyAny>, x: f64, @@ -128,7 +128,7 @@ pub fn third_partial_derivative( /// second partial derivative w.r.t. variables i and k /// second partial derivative w.r.t. variables j and k /// third partial derivative -#[allow(clippy::type_complexity)] +#[expect(clippy::type_complexity)] pub fn third_partial_derivative_vec( f: &Bound<'_, PyAny>, x: Vec, From c8a6a3fde01b8e754c058fb23867c184be2e3a48 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Thu, 5 Sep 2024 20:15:19 +0200 Subject: [PATCH 2/2] add rust version requirement --- Cargo.toml | 1 + src/dual_size.rs | 309 ----------------------------------------------- 2 files changed, 1 insertion(+), 309 deletions(-) delete mode 100644 src/dual_size.rs diff --git a/Cargo.toml b/Cargo.toml index 5b722f0..da5b3d0 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -3,6 +3,7 @@ name = "num-dual" version = "0.9.1" authors = ["Gernot Bauer ", "Philipp Rehner "] +rust-version = "1.81" edition = "2021" readme = "README.md" license = "MIT OR Apache-2.0" diff --git a/src/dual_size.rs b/src/dual_size.rs deleted file mode 100644 index 2fe0bee..0000000 --- a/src/dual_size.rs +++ /dev/null @@ -1,309 +0,0 @@ -use crate::{DualNum, DualNumFloat, OneZero, Scale}; -use nalgebra::allocator::Allocator; -use nalgebra::*; -use num_traits::{Float, Inv, One}; -use std::convert::Infallible; -use std::fmt; -use std::iter::{Product, Sum}; -use std::marker::PhantomData; -use std::ops::{ - Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Rem, RemAssign, Sub, SubAssign, -}; - -/// A dual number for the calculations of gradients or Jacobians. -#[derive(PartialEq, Eq, Clone, Debug)] -pub struct DualVec, F, N: Dim> -where - DefaultAllocator: Allocator, -{ - /// Real part of the dual number - pub re: T, - /// Derivative part of the dual number - pub eps: OVector, - f: PhantomData, -} - -impl + Copy, F: Copy, const N: usize> Copy for DualVec> {} - -pub type DualVec32 = DualVec; -pub type DualVec64 = DualVec; -pub type DualSVec32 = DualVec>; -pub type DualSVec64 = DualVec>; -pub type DualDVec32 = DualVec; -pub type DualDVec64 = DualVec; -pub type Dual = DualVec; -pub type Dual32 = Dual; -pub type Dual64 = Dual; - -impl, F, N: Dim> DualVec -where - DefaultAllocator: Allocator, -{ - /// Create a new dual number from its fields. - #[inline] - pub fn new(re: T, eps: OVector) -> Self { - Self { - re, - eps, - f: PhantomData, - } - } -} - -impl, F, N: Dim> OneZero for DualVec -where - DefaultAllocator: Allocator, -{ - fn zero(&self) -> Self { - let (n, _) = self.eps.shape_generic(); - Self::new( - self.re.zero(), - OVector::from_element_generic(n, U1, self.re.zero()), - ) - } - - fn one(&self) -> Self { - let (n, _) = self.eps.shape_generic(); - Self::new( - self.re.one(), - OVector::from_element_generic(n, U1, self.re.zero()), - ) - } -} - -impl, F> Dual { - /// Create a new scalar dual number from its fields. - #[inline] - pub fn new_scalar(re: T, eps: T) -> Self { - Self::new(re, SVector::from_element(eps)) - } -} - -// impl + Zero, F, N: Dim> DualVec -// where -// DefaultAllocator: Allocator, -// { -// /// Create a new dual number from the real part. -// #[inline] -// pub fn from_re(re: T) -> Self { -// Self::new(re, OVector::zeros_generic()) -// } -// } - -impl + One, F> Dual { - /// Set the derivative part to 1. - /// ``` - /// # use num_dual::{Dual64, DualNum}; - /// let x = Dual64::from_re(5.0).derivative().powi(2); - /// assert_eq!(x.re, 25.0); - /// assert_eq!(x.eps.unwrap(), 10.0); - /// ``` - #[inline] - pub fn derivative(x: T) -> Self { - let one = x.one(); - Self::new(x, SVector::from_element(one)) - } -} - -/// Calculate the first derivative of a scalar function. -/// ``` -/// # use num_dual::{first_derivative, DualNum}; -/// let (f, df) = first_derivative(|x| x.powi(2), 5.0); -/// assert_eq!(f, 25.0); -/// assert_eq!(df, 10.0); -/// ``` -pub fn first_derivative, F>(g: G, x: T) -> (T, T) -where - G: FnOnce(Dual) -> Dual, -{ - try_first_derivative(|x| Ok::<_, Infallible>(g(x)), x).unwrap() -} - -/// Variant of [first_derivative] for fallible functions. -pub fn try_first_derivative, F, E>(g: G, x: T) -> Result<(T, T), E> -where - G: FnOnce(Dual) -> Result, E>, -{ - let one = x.one(); - let x = Dual::new_scalar(x, one); - g(x).map(|r| (r.re, r.eps.data.0[0][0].clone())) -} - -/// Calculate the gradient of a scalar function -/// ``` -/// # use approx::assert_relative_eq; -/// # use num_dual::{gradient, DualNum, DualSVec64}; -/// # use nalgebra::SVector; -/// let v = SVector::from([4.0, 3.0]); -/// let fun = |v: SVector, 2>| (v[0].powi(2) + v[1].powi(2)).sqrt(); -/// let (f, g) = gradient(fun, v); -/// assert_eq!(f, 5.0); -/// assert_relative_eq!(g[0], 0.8); -/// assert_relative_eq!(g[1], 0.6); -/// ``` -/// -/// The variable vector can also be dynamically sized -/// ``` -/// # use approx::assert_relative_eq; -/// # use num_dual::{gradient, DualNum, DualDVec64}; -/// # use nalgebra::DVector; -/// let v = DVector::repeat(4, 2.0); -/// let fun = |v: DVector| v.iter().map(|v| v * v).sum::().sqrt(); -/// let (f, g) = gradient(fun, v); -/// assert_eq!(f, 4.0); -/// assert_relative_eq!(g[0], 0.5); -/// assert_relative_eq!(g[1], 0.5); -/// assert_relative_eq!(g[2], 0.5); -/// assert_relative_eq!(g[3], 0.5); -/// ``` -pub fn gradient, F: DualNumFloat, N: Dim>( - g: G, - x: OVector, -) -> (T, OVector) -where - G: FnOnce(OVector, N>) -> DualVec, - DefaultAllocator: Allocator + Allocator, N>, -{ - try_gradient(|x| Ok::<_, Infallible>(g(x)), x).unwrap() -} - -/// Variant of [gradient] for fallible functions. -pub fn try_gradient, F: DualNumFloat, E, N: Dim>( - g: G, - x: OVector, -) -> Result<(T, OVector), E> -where - G: FnOnce(OVector, N>) -> Result, E>, - DefaultAllocator: Allocator + Allocator, N>, -{ - let (n, _) = x.shape_generic(); - let x = x.map_with_location(|i, _, x| { - let mut eps = OVector::from_element_generic(n, U1, x.zero()); - eps[i] = x.one(); - DualVec::new(x, eps) - }); - g(x).map(|res| (res.re, res.eps)) -} - -/// Calculate the Jacobian of a vector function. -/// ``` -/// # use num_dual::{jacobian, DualSVec64, DualNum}; -/// # use nalgebra::SVector; -/// let xy = SVector::from([5.0, 3.0, 2.0]); -/// let fun = |xy: SVector, 3>| SVector::from([ -/// xy[0] * xy[1].powi(3) * xy[2], -/// xy[0].powi(2) * xy[1] * xy[2].powi(2) -/// ]); -/// let (f, jac) = jacobian(fun, xy); -/// assert_eq!(f[0], 270.0); // xy³z -/// assert_eq!(f[1], 300.0); // x²yz² -/// assert_eq!(jac[(0,0)], 54.0); // y³z -/// assert_eq!(jac[(0,1)], 270.0); // 3xy²z -/// assert_eq!(jac[(0,2)], 135.0); // xy³ -/// assert_eq!(jac[(1,0)], 120.0); // 2xyz² -/// assert_eq!(jac[(1,1)], 100.0); // x²z² -/// assert_eq!(jac[(1,2)], 300.0); // 2x²yz -/// ``` -pub fn jacobian, F: DualNumFloat, M: Dim, N: Dim>( - g: G, - x: OVector, -) -> (OVector, OMatrix) -where - G: FnOnce(OVector, N>) -> OVector, M>, - DefaultAllocator: Allocator, M> - + Allocator - + Allocator - + Allocator - + Allocator, N> - + Allocator, N> - + Allocator, M>, -{ - try_jacobian(|x| Ok::<_, Infallible>(g(x)), x).unwrap() -} - -/// Variant of [jacobian] for fallible functions. -#[allow(clippy::type_complexity)] -pub fn try_jacobian, F: DualNumFloat, E, M: Dim, N: Dim>( - g: G, - x: OVector, -) -> Result<(OVector, OMatrix), E> -where - G: FnOnce(OVector, N>) -> Result, M>, E>, - DefaultAllocator: Allocator, M> - + Allocator - + Allocator - + Allocator - + Allocator, N> - + Allocator, N> - + Allocator, M>, -{ - let (n, _) = x.shape_generic(); - let x = x.map_with_location(|i, _, x| { - let mut eps = OVector::from_element_generic(n, U1, x.zero()); - eps[i] = x.one(); - DualVec::new(x, eps) - }); - g(x).map(|res| { - let eps = OMatrix::from_rows(res.map(|res| res.eps.transpose()).as_slice()); - (res.map(|r| r.re), eps) - }) -} - -/* chain rule */ -impl, F: Float, N: Dim> DualVec -where - DefaultAllocator: Allocator, -{ - #[inline] - fn chain_rule(&self, f0: T, f1: T) -> Self { - Self::new(f0, self.eps.clone() * f1) - } -} - -/* product rule */ -impl<'a, 'b, T: DualNum, F: Float, N: Dim> Mul<&'a DualVec> for &'b DualVec -where - DefaultAllocator: Allocator, -{ - type Output = DualVec; - #[inline] - fn mul(self, other: &DualVec) -> Self::Output { - DualVec::new( - self.re.clone() * other.re.clone(), - &self.eps * other.re.clone() + &other.eps * self.re.clone(), - ) - } -} - -/* quotient rule */ -impl<'a, 'b, T: DualNum, F: Float, N: Dim> Div<&'a DualVec> for &'b DualVec -where - DefaultAllocator: Allocator, -{ - type Output = DualVec; - #[inline] - fn div(self, other: &DualVec) -> DualVec { - let inv = other.re.recip(); - DualVec::new( - self.re.clone() * inv.clone(), - (self.eps.clone() * other.re.clone() - other.eps.clone() * self.re.clone()) - * inv.clone() - * inv, - ) - } -} - -/* string conversions */ -impl, F, N: Dim> fmt::Display for DualVec -where - DefaultAllocator: Allocator, -{ - fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { - todo!() - // write!(f, "{}", self.re)?; - // self.eps.fmt(f, "ε") - } -} - -impl_first_derivatives2!(DualVec, [eps], [D]); -impl_dual2!(DualVec, [eps], [D]);