Add nalgebra compatibility for DualVec #59
Conversation
|
Here's an atan2 on hyper-duals at least https://github.com/oberbichler/HyperJet/blob/main/include/hyperjet/hyperjet.h#L1412 |
|
Thank you for the very nice contribution! This should also enable us to switch from the custom implementation of eigenvalues to nalgebra. I will finish work on #58 which was basically only waiting for confirmation that downstream crates do not suffer from any performance loss. Afterwards, I would ask you to rebase this branch and I hope to have a closer look at this soon. |
|
@cormacrelf, this is fantastic! I've been working on implementing this exact integration but yours is far superior to mine so I will cease development of my work and await this merge to achieve the functionality that I desire 😁 I feel your pain around the SimdValue issues as I spent a long time banging my head against that wall before resigning myself to the ComplexField way! @prehner, do you envisage this crate ever extending beyond reals to complex? Or is this something that is beyond the scope of the package? (Since that would entail some form of quaternion algebra) I ask since the nalgebra integration here assumes, as it rightly should, that F is bound to Float such that the ComplexField has no imaginary part. If I can help in anyway with regards to review or contribution then please let me know! P.S. For the atan2 of duals, an implementation might look something like this: fn atan2(self, other: Self) -> Self {
let re = self.re.atan2(other.re);
let dx = (self.re * other.eps.into()) / (self.re.powi(2) + other.re.powi(2));
let dy = (other.re * self.eps.into()) / (self.re.powi(2) + other.re.powi(2));
DualVec {
re: re,
eps: StaticVec::from(dy-dx),
f: PhantomData, // Is this aspect correct?
}
} |
|
With regards to extending to complex numbers: At the moment, I think I don't see where exactly the problem comes from, but that is likely due to me also giving up on the complexity of the nalgebra trait system. In theory you could have a dual number with complex numbers as elements (that would require the complex type to implement DualNum for most of the operations) or have a complex number with dual numbers as elements (Then the implementation of the complex number data type has to be generic enough to allow that). Mathematically that is the same. With regards to atan2: Your suggested implementation looks good. The PhantomData can be hidden by simply calling the |
|
Yes. And for example the test I added in this PR uses a UnitQuaternion. Quaternions are complex numbers with three imaginary parts, but like the 2-dimensional complex numbers, all four parts are represented as reals. The imaginary bits are just coefficients . This crate's types are perfectly capable of emulating a real number for nalgebra by implementing the RealField or SimdRealField trait. The confusion around this probably comes from simba/nalgebra's traits not having enough docs, ie RealField requiring ComplexField as well without really explaining why. I note that the implementation of ComplexField for this crate follows the one for f32 and f64. There simply isn't a complex part. It returns zero when you ask and otherwise behaves as if there is no complex part (so norm1() == abs()). AFAICT The trait ComplexField does not map to a mathematical trait like "field" or "ring", it is merely a way to make stuff like "multiplying a It's certainly possible that nobody has wanted to implement these traits enough to warrant documenting them. Certainly the issue I linked makes it seem that way, as you would think someone would have found it before. Since I have wrapped my head around it, I should probably send some docs PRs to simba. |
|
The problems that I imagined might occur are around the interaction between the imaginary and dual units when performing the algebra as a result of a complex dual being a quaternion: For example on page 3, division is not defined for complex duals and complex duals have five distinct types of conjugation. As a result, the reciprocal of a complex dual with no real part is undefined and there are five different types of possible norm with their own nuances. Perhaps this is achievable with a custom-defined crate-bound Complex type but this discussion is beyond the scope of this PR so I'll drop it! |
|
Division is only a problem if the real part is zero, which comes at no surprise. Some methods, like the norm, might have different meanings/results for a @cormacrelf Can you rebase this PR on the current master? |
|
Yep, just took a lot of wrangling and I've been busy with other things. Everything needed truly awful amounts of |
fixup maybeuninit, zero
…:pi() etc with it
Update commentary on min/max/clamp Better panic for copysign and atan2 implement complex norm1 using abs
eps[0] to eps.unwrap()
|
Looking forward to having this functionality merged, thanks again both! |
There was a problem hiding this comment.
Thank you for the patience. The rebase and the introduction of the Derivative struct goes indeed quite against the SIMD implementation. After all it was done to allow for dynamically sized dual numbers which are probably not optimal for SIMD. However, as you write this is mostly about unlocking nalgebra functionalities that are hidden behind these traits and not actually optimizing dual numbers with SIMD.
Going through this PR really reinforced my suspicion that nalgebra has design flaws when it comes to trait bounds and the usage of data types aside from f32 and f64. But how are the chances that they tighten it up when that is not the main focus of the crate?
|
|
||
| mod nalgebra_api { | ||
| use nalgebra::{Point2, Point3, UnitQuaternion, Vector2, Vector3}; | ||
| use num_dual::*; | ||
| use num_traits::Zero; | ||
| use std::f32::consts::*; | ||
|
|
||
| fn unit_circle(t: Dual32) -> Point2<Dual32> { | ||
| let x_dir = |x: Dual32| Vector2::new(x, Dual32::zero()); | ||
| let y_dir = |y: Dual32| Vector2::new(Dual32::zero(), y); | ||
| let theta = t * TAU; | ||
| Point2::from(x_dir(theta.cos()) + y_dir(theta.sin())) | ||
| } | ||
|
|
||
| // This is testing that you can type-check code that whacks DualVec in | ||
| // nalgebra structures and tries to use them. | ||
| #[test] | ||
| fn use_nalgebra_2d() { | ||
| // 1 radian around the circle | ||
| let t = Dual32::from_re(0.25).derivative(); | ||
| let point = unit_circle(t); | ||
| let real = point.map(|x| x.re); | ||
| let grad = point.map(|x| x.eps.unwrap()); | ||
| println!("point: {}", point.coords); | ||
| approx::assert_relative_eq!(real, Point2::new(0., 1.), epsilon = 1e-3); | ||
| approx::assert_relative_eq!(grad, Point2::new(-TAU, 0.), epsilon = 1e-3); | ||
| } | ||
|
|
||
| #[test] | ||
| fn use_nalgebra_3d() { | ||
| // First one does nothing to the gradient. Still got no y or z direction. | ||
| let rot1 = UnitQuaternion::from_axis_angle(&Vector3::x_axis(), FRAC_PI_8); | ||
| // Second one goes pi/8 (22.5deg) about the y axis, which should shift some of the steepness | ||
| // in the x direction into the z direction, but not much. | ||
| let rot2 = UnitQuaternion::from_axis_angle(&Vector3::y_axis(), FRAC_PI_8); | ||
| let rotation = (rot2 * rot1).cast::<Dual32>(); | ||
| let lifted_3d_circle = |t: Dual32| { | ||
| let xy = unit_circle(t); // [0, 1] with derivatives [-1, 0] | ||
| Point3::new(xy.x, xy.y, Dual32::zero()) | ||
| }; | ||
| let function = |t: Dual32| rotation * lifted_3d_circle(t); | ||
| let point = function(Dual32::from_re(0.25).derivative()); | ||
| let real = point.map(|x| x.re); | ||
| let grad = point.map(|x| x.eps.unwrap()); | ||
| println!("rotated point: {}", point.coords); | ||
| approx::assert_relative_eq!(real.coords, real.coords.normalize(), epsilon = 1e-3); | ||
| approx::assert_relative_eq!(real, Point3::new(0.146, 0.924, 0.354), epsilon = 1e-3); | ||
| approx::assert_relative_eq!(grad, Point3::new(-5.805, 0., 2.404), epsilon = 1e-3); | ||
| } | ||
| } |
There was a problem hiding this comment.
This file is auto generated by the jupyter notebook in the repo. Please move to a separate file.
| Self { | ||
| re, | ||
| eps, | ||
| f: PhantomData, | ||
| } |
There was a problem hiding this comment.
For readability, I suggest using the new constructor in these places instead.
| @@ -4,7 +4,9 @@ use nalgebra::constraint::{SameNumberOfRows, ShapeConstraint}; | |||
| use nalgebra::*; | |||
There was a problem hiding this comment.
clippy has some suggestions in this file that you might want to consider.
Does what it says on the tin. Partially addresses #55
0.35355338 + 2.404471εThe Simd-enabled part turned out to be useless for now, see dimforge/simba#44. It was pretty easy to implement SimdValue and thus we are able to split a
DualVec<simba::simd::f32x4, f32, 1>into fourDualVec<f32, f32, 1>, and join them back together. So we can get in and out of SIMD-land.Alas, we are not able to implement SimdRealField/SimdComplexField directly and thus use T's vectorized operations. So we are stuck with implementing the non-vectorized
RealField/ComplexFieldfor now. If that problem is fixed, we'll be able to swap out that with the Simd version of the trait, using basically find/replacefnwithfn simd_or thereabouts.