search

itt-ustutt / num-dual

Generalized (hyper-) dual numbers in rust
Other
53 stars 6 forks source link

Add nalgebra compatibility for DualVec #59

Closed cormacrelf closed 1 year ago

cormacrelf commented 1 year ago

Does what it says on the tin. Partially addresses #55

The Simd-enabled part turned out to be useless for now, see https://github.com/dimforge/simba/issues/44. It was pretty easy to implement SimdValue and thus we are able to split a DualVec<simba::simd::f32x4, f32, 1> into four DualVec<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/ComplexField for now. If that problem is fixed, we'll be able to swap out that with the Simd version of the trait, using basically find/replace fn with fn simd_ or thereabouts.

cormacrelf commented 1 year ago

Here's an atan2 on hyper-duals at least

https://github.com/oberbichler/HyperJet/blob/main/include/hyperjet/hyperjet.h#L1412

cormacrelf commented 1 year ago

Atan2 on duals https://github.com/LLNL/serac/blob/develop/src/serac/numerics/functional/dual.hpp#L328

prehner commented 1 year ago

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.

tsmith023 commented 1 year ago

@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 :grin: 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?
    }
}
prehner commented 1 year ago

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 new constructor.

cormacrelf commented 1 year ago

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 Complex<f32> by an f32" a bit faster since you don't need to upgrade the f32 first, it already knows how to present itself as a complex number.

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.

tsmith023 commented 1 year ago

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: $$\huge z+t\varepsilon=x+yi+p\varepsilon+qi\varepsilon$$ of the form $\large a+bi+cj+dk$ with $\large i^2=-1$, $\large j^2=0$, and $\large k^2=0$.

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!

prehner commented 1 year ago

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 Dual<Complex<f64>> and a Complex<Dual<f64>> but I don't see that as a big problem. After all, in practice that would depend on what the actual application is. Now that I think about it, we actually use Complex<T> with T: DualNum<f64> in some part of our code. Both the real and imaginary parts of the complex number then depend on some input variable and the dual parts of each contain the derivatives with respect to that variable.

@cormacrelf Can you rebase this PR on the current master?

cormacrelf commented 1 year ago

Yep, just took a lot of wrangling and I've been busy with other things. Everything needed truly awful amounts of DefaultAllocator trait bounds (more than you would think!) and the move to the Derivative with an Option in it was painful. But it worked out.

tsmith023 commented 1 year ago

Looking forward to having this functionality merged, thanks again both!