JuliaManifolds / Manifolds.jl

Manifolds.jl provides a library of manifolds aiming for an easy-to-use and fast implementation.
https://juliamanifolds.github.io/Manifolds.jl
MIT License
368 stars 53 forks source link

Minor bug with `rand` at `Identity` #665

Closed olivierverdier closed 10 months ago

olivierverdier commented 10 months ago

I'm trying the very last version of the Manifolds ecosystem (love it so far!), and the following code that worked before doesn't anymore:

G = SpecialOrthogonal(3)
rand(G; vector_at=Identity(G))

It's not a big deal, as rand(G; vector_at=identity_element(G)) does work, but perhaps it could hide a potentially more harmful bug? If not, feel free to close the issue right away.

kellertuer commented 10 months ago

We changed random tangent vector generation slightly, since for non-array-represetnations (struct) we had to fix the tangent vector allocation. I am quite sure it does not hide a potentially more harmful bug, we just missed to transform the Identity(G) type (which does not allocate a full identity matrix here but is a “memory-efficient” struct) to something where we can allocate the array correctly for the tangent vector (Lie algebra element) here.

But you are right, reintroducing that would surely be nice.

mateuszbaran commented 10 months ago

Thanks for reporting it. I've made a PR that fixes this issue for some groups (I don't really see a good way to solve it without doing group-specific additions).

olivierverdier commented 10 months ago

I have a suggestion regarding that: what about making rand(TangentSpace(G, Identity(G))) work instead? (and deprecate rand(G; vector_at=...)?).

The advantage: the error message for a missing such rand methods will be easier on the user. Either with Method error no method matching .... rand(... Identity), or you could use a fallback method such as rand(::Fiber{<:Any, <:Any, <:Any, <:Identity}) = error("Not defined, useidentity_elementinstead").. Does that make sense?

mateuszbaran commented 10 months ago

I'm not sure about deprecating rand(G; vector_at=...) but rand(TangentSpace(G, Identity(G))) looks like something that should work :+1:

kellertuer commented 10 months ago

Yes I think we introduced the TangentSpace one with the last rewrite as well. But I would also prefer to not deprecate the vector_at format, since not everyone wants to think of generating a tangent space first (and the tangent space one does fall back to the verctor_at case for now).

olivierverdier commented 10 months ago

Then it would be more useful to the user to have a dedicated function, like random_vector_at(M, x) = rand(TangentSpace(M,x)).

Again, the main advantage is that you can implement a general fallback, or have much better error messages, when rand is not implemented for Identity, across the whole library: random_vector(::Any, ::Identity) = error(...).

I also think it's easier to use and remember than the vector_at keyword...

kellertuer commented 10 months ago

We discussed that quite a while back (maybe a few years) where we indeed hat rand/rand! and random_vector(I think) was the name. The problem is that this does not fit with the naming scheme of random functions in the (global) Julia ecosystem. Similarly (if we would adapt our notation) something like rand(M, p, X).

In the current form, the kwargs...we accept and pass on (for example when generating a random matrix somewhere) are consistent and fit the rand scheme of Julia.

We also abandoned throwing our own error messages, because the MethodErrors usually propose more alternatives that do exist

olivierverdier commented 10 months ago

Mmh, I cannot see that rand uses keyword arguments neither in Base nor Distributions, but I may be wrong... On the other hand I see that there is randexp, randstring, randsubseq, randperm, etc., so randvector could fit pretty well.

As to the errors, I agree. As I said above, you would get a nice Method error no method matching .... randvector(... Identity), which is better than the current obscure error messages when using Identity instead of identity_element.

I'll leave it to you to decide what is best for the users. ☺️

kellertuer commented 10 months ago

we had rand_vector a few years back. Besides such a change would be a (heavily breaking) change – we specifically decided for one rand (rand!) command back then.

The error above is still not one we throw but due to the decorator pattern we employ (in order to implement several metrics for one manifold for example). On the other hand Identity is a very specific type, that might not yet be fully tested, because until now there seems to have been not much interest in the Lie group part so it exists but was not used much, neither by the developers nor by other users. Then bugs might exist.