Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Power Manifolds usage and bugs #497

Closed
Affie opened this issue Jul 4, 2022 · 23 comments
Closed

Power Manifolds usage and bugs #497

Affie opened this issue Jul 4, 2022 · 23 comments

Comments

@Affie
Copy link
Contributor

Affie commented Jul 4, 2022

Hello,
As Mateusz sugested in #495 (comment)
I experimented with using power manifolds to represent points for optimization.

From the documentation, I figured the NestedReplacingPowerRepresentation with an ArrayPartition of <:StaticArrays looked the most promising. So I tried it and ran into errors, so I tried NestedPowerRepresentation with ArrayPartition and again errors.
I completely understand Manifolds.jl is still under development so I tried to help by writing a simple script for testing which combination to try.

I ended up using "ProductRepr SA NestedReplacingPowerRepresentation" and I got it to work.

  1. What is your recommendation of what PowerManifold to use? So we can focus on implementing it and fixing bugs as they come up.
  2. What is your recommendation between ProductRepr and ArrayPartition?
    I think you already answered it here (so linking for the next person): There are almost no drawbacks to ArrayPartition vs ProductRepr
  3. How do I do group operations such as, identity_element,vee,hat
  4. I think a PowerManifold dispatch such as log(M, Identity(M), p) where all points are calculated from the identity element can be quite handy, I currently create a vector of identity elements and use that.

Thanks again for the support. I tried to fix the errors myself, but couldn't quite figure out the trait macros.

Here is the script in case you find it handy:

using Manifolds
using RecursiveArrayTools: ArrayPartition
using StaticArrays
using Test
using PreallocationTools

G = SpecialEuclidean(2)

# the different options
options = []
#ProductRepr
o = (name = "ProductRepr NestedPowerRepresentation",
     p1 = ProductRepr([0.0,0.0], [1.0 0.0; 0.0 1.0]),
     p2 = ProductRepr([0.0,0.0], [0.0 -1.0; 1.0 0.0]),
     M = PowerManifold(G, NestedPowerRepresentation(), 2));
push!(options, o);

#ArrayPartition
o = (name = "ArrayPartition NestedPowerRepresentation",
     p1 = ArrayPartition([0.0,0.0], [1.0 0.0; 0.0 1.0]),
     p2 = ArrayPartition([0.0,0.0], [0.0 -1.0; 1.0 0.0]),
     M = PowerManifold(G, NestedPowerRepresentation(), 2));
push!(options, o);

#ProductRepr with NestedReplacingPowerRepresentation
o = (name = "ProductRepr NestedReplacingPowerRepresentation",
     p1 = ProductRepr([0.0,0.0], [1.0 0.0; 0.0 1.0]), 
     p2 = ProductRepr([0.0,0.0], [0.0 -1.0; 1.0 0.0]),
     M = PowerManifold(G, NestedReplacingPowerRepresentation(), 2));
push!(options, o);

o = (name = "ProductRepr SA NestedReplacingPowerRepresentation",
     p1 = ProductRepr(SA[0.0,0.0], SA[1.0 0.0; 0.0 1.0]),
     p2 = ProductRepr(SA[0.0,0.0], SA[0.0 -1.0; 1.0 0.0]),
     M = PowerManifold(G, NestedReplacingPowerRepresentation(), 2));
push!(options, o);

#ArrayPartition
o = (name = "ArrayPartition NestedReplacingPowerRepresentation",
     p1 = ArrayPartition([0.0,0.0], [1.0 0.0; 0.0 1.0]),
     p2 = ArrayPartition([0.0,0.0], [0.0 -1.0; 1.0 0.0]),
     M = PowerManifold(G, NestedReplacingPowerRepresentation(), 2));
push!(options, o);

o = (name = "ArrayPartition SA NestedReplacingPowerRepresentation",
     p1 = ArrayPartition(SA[0.0,0.0], SA[1.0 0.0; 0.0 1.0]),
     p2 = ArrayPartition(SA[0.0,0.0], SA[0.0 -1.0; 1.0 0.0]),
     M = PowerManifold(G, NestedReplacingPowerRepresentation(), 2));
push!(options, o);

#ArrayPartition
o = (name = "Vector NestedReplacingPowerRepresentation",
     p1 = [0.0,0.0],
     p2 = [1.0,1.0],
     M = PowerManifold(TranslationGroup(2), NestedReplacingPowerRepresentation(), 2));
push!(options, o);

o = (name = "SVector NestedReplacingPowerRepresentation",
     p1 = SA[0.0,0.0],
     p2 = SA[1.0,1.0],
     M = PowerManifold(TranslationGroup(2), NestedReplacingPowerRepresentation(), 2));
push!(options, o);

## Try to find one that fully Works


for (i,o) in enumerate(options)
  @info "option $i: $(o.name)"
  try
     p = [o.p1,o.p1]
     q = [o.p1,o.p2]
     M = o.M
     log(M, p, q)

     X = allocate(p)

     log!(M, X, p, q)
     
     X = [allocate(o.p1),allocate(o.p1)]
     log!(M, X, p, q)

     exp(M, p, X)

     exp!(M, q, p, X)

     Xc = get_coordinates(M, p, X, DefaultBasis())

     get_coordinates!(M, Xc, p, X, DefaultBasis())

     get_vector(M, p, Xc, DefaultBasis())

     get_vector!(M, X, p, Xc, DefaultBasis())

     # @info "Success of Manifols functions, checking dualcache"
     # dualcache(X)
     @info "Success"
     
     println()
  catch ex
     @error ex
     println()
  end
end
[ Info: option 1: ProductRepr NestedPowerRepresentation
┌ Error: MethodError(view, (ProductRepr{Tuple{Vector{Float64}, Matrix{Float64}}}(([6.90992505373676e-310, 6.9099250553581e-310], [0.0 6.9099139117185e-310; 6.90991392036585e-310 0.0])), Colon()), 0x0000000000007d1e)
└ 

[ Info: option 2: ArrayPartition NestedPowerRepresentation
┌ Error: ErrorException("type SubArray has no field parts")
└ 

[ Info: option 3: ProductRepr NestedReplacingPowerRepresentation
[ Info: Success

[ Info: option 4: ProductRepr SA NestedReplacingPowerRepresentation
[ Info: Success

[ Info: option 5: ArrayPartition NestedReplacingPowerRepresentation
[ Info: Success

[ Info: option 6: ArrayPartition SA NestedReplacingPowerRepresentation
┌ Error: MethodError(convert, (ArrayPartition{Float64, Tuple{SVector{2, Float64}, SMatrix{2, 2, Float64, 4}}}, ArrayPartition{Float64, Tuple{SVector{2, Float64}, MMatrix{2, 2, Float64, 4}}}(([0.0, 0.0], [0.0 -0.0; 0.0 0.0]))), 0x0000000000007d1e)
└ 

[ Info: option 7: Vector NestedReplacingPowerRepresentation
[ Info: Success

[ Info: option 8: SVector NestedReplacingPowerRepresentation
[ Info: Success
@kellertuer
Copy link
Member

kellertuer commented Jul 4, 2022

Hm, in my mind neither ProductRepr nor ArrayPartition are meant for Power manifolds, they are for product manifolds (as the first might indicate ;) ).

Power manifolds are far more homogeneous in the sense that on M^3 we know that we have “three times the same thing”. So there is three things how to work with PowerManifold (see https://juliamanifolds.github.io/Manifolds.jl/latest/manifolds/power.html) - I try to illustrate them on M = PowerManifold(Sphere(2),3)

  • ArrayPowerRepresentation, that is you have one array containing the values, the first array dimension is the manifold (for the sphere, might be two dimensions for matrices) the next (second or third) is the dimension for the power. In other words
p = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 -1.0]

is three columns of points on the manifold and p[M,3] returns the last column [0.0;0.0;-1.0]

My disadvantage in this is for images of matrix-valued data (4D arrays) this is a little clumsy (though I only worked with this in Matlab and it was my main reason to switch)

bith nested ones are similar but the same point as above would look like

M2 = PowerManifold(Sphere(2), NestedPowerRepresentation(), 3)
p2 = [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0,0.0, -1.0]]

Again of course p2[M2,3] works the same as above (one reason we introduced this array indexing thingy)

so here we have an array (of length 3, since we have length 3 of power manifold) where each entry is an array itself of a point on S2. So it is nested arrays. The replacing one is just a variant of this handling the internal arrays in the nested setting differently (for technical reasons that is sometimes nicer).

  1. How do I do group operations such as, identity_element,vee,hat

A power manifold does not have a group structure per se, so that again depends heavily on what you use as M. What we could do is for the case that M is a group then implicitly the power manifold is a group as well (kind of the power group), the new system (ManifoldsBase.jl 0.13) allows for this automatic interplay for sure, but it is not implemented (yet).

  1. I think a PowerManifold dispatch such as log(M, Identity(M), p) where all points are calculated from the identity element can be quite handy, I currently create a vector of identity elements and use that.

Similar to 3, we first need an identity on a power manifold (given that it is a group, for example by element wise group action which would be the default) then the identity would exist automatically and the log would work element wise (on the identities as well, sure)

PS: I wouldn't‘t say that the general Manifolds interface is that much under development anymore, the new (0.13.0 ManifoldsBase) one is quite robust and nice I hope. What is not yet covered that much (because we do not have many people using it) is the Lie group part, especially when its about the interplay between product, power, metric and Lie. This is also quite tricky as quite some discussions show, since there is usually not one single way to decide for a certain default.

@Affie
Copy link
Contributor Author

Affie commented Jul 4, 2022

Hm, in my mind neither ProductRepr nor ArrayPartition are meant for Power manifolds, they are for product manifolds (as the first might indicate ;) ).

Sorry, I should have been clearer. I'm currently only testing with Points (TranslationGroup) and Poses (SE)

The problem can have a lot of points and Poses. So I figured I'd start with

M1 = PowerManifold(TranslationGroup(N), NestedReplacingPowerRepresentation(), number_of_points_on_1)
M2 = PowerManifold(SpecialEuclidean(N), NestedReplacingPowerRepresentation(), number_of_points_on_2)
M = ProductManifold(M1, M2)

Then points on M can be ProductRepr or ArrayPartition
The points on SpecialEuclidean in this case can also be ProductRepr or ArrayPartition or just a Matrix

@Affie
Copy link
Contributor Author

Affie commented Jul 4, 2022

A power manifold does not have a group structure per se, so that again depends heavily on what you use as M. What we could do is for the case that M is a group then implicitly the power manifold is a group as well (kind of the power group), the new system (ManifoldsBase.jl 0.13) allows for this automatic interplay for sure, but it is not implemented (yet).

I think that is what I'm after, thanks

@kellertuer
Copy link
Member

M1 = PowerManifold(TranslationGroup(N), NestedReplacingPowerRepresentation(), number_of_points_on_1)
M2 = PowerManifold(SpecialEuclidean(N), NestedReplacingPowerRepresentation(), number_of_points_on_2)
M = ProductManifold(M1, M2)

Ah this is quite a complicated structure then, sure, the first is just an array of arrays, since translation group is just Rn, the second is an array of ProductRepr? And those two should be part of a ProductRepr then.

The main problem here (again) is that the power manifold does not “look into” its (element) manifold to inherit its poperties (in principle it could, since 0.13.0 automatically become a PowerLieGroup for example, but for now it does not). The thing here is that for any of these interplays (of power with something or Product with something - like metric, like Lie like any other decorated feature) this is only doable in a reasonable way since 0.13.0 and I was happy to finish that in general, so these interplays are a next step.

  • when are they reasonable?
  • does Power(Lie(M)) always have a Lie structure? Is it the same as Lie(Power(M))?
  • ...

To me this is not so clear but if this is clear, we can surely provide that

One general remark why I said ProductRepr is the wrong thing on Power manifolds:

I now know what confused me. Here

G = SpecialEuclidean(2)
p1 = ProductRepr([0.0,0.0], [1.0 0.0; 0.0 1.0])
p2 = ProductRepr([0.0,0.0], [0.0 -1.0; 1.0 0.0])
M = PowerManifold(G, NestedReplacingPowerRepresentation(), 2)

It looked like p1 and p2 where points on M, but it was a long day and you actually do p = [p1, p2] later.
Then I am not sure why ArrayPartition is failing, but nice that the rest works fine (I am also more a PowerManifolds user, Mateusz is the Product Expert for sure)

@mateuszbaran
Copy link
Member

  1. What is your recommendation of what PowerManifold to use? So we can focus on implementing it and fixing bugs as they come up.

Yes, we definitely need more documentation about which types of power manifold are suitable for which cases. "ArrayPartition SA NestedReplacingPowerRepresentation" should work fine in this case, I will look into why it fails.

Ronny has explained quite nicely how it works, and I also think that the interface works well here and the issues are that not everything is implemented and documented, and such complex interactions between power, product and Lie need to be ironed out -- but that definitely looks doable to me. I will look into implementing some missing things tomorrow.

@mateuszbaran
Copy link
Member

#498 and JuliaManifolds/ManifoldsBase.jl#113 should make all these options working. ArrayPartition SA NestedReplacingPowerRepresentation should be overall the best but it requires some optimization, I will continue with optimizing what you need.

@Affie
Copy link
Contributor Author

Affie commented Jul 5, 2022

Thanks a lot. We will work towards 'ArrayPartition SA NestedReplacingPowerRepresentation' and I'll open separate issues as they come up.
Maybe one day I'll be able to open PR to help fix issues, but for now, I'll help you make Manifolds better by using it and logging what I find.

@mateuszbaran
Copy link
Member

You're welcome 🙂 . Your reports are really helpful (they tell me which parts need improvements) and I appreciate that you are willing to deal with these issues.

mateuszbaran added a commit that referenced this issue Jul 7, 2022
* Fixing some issues from #497

* some non-mutating methods for Rotations(2)

* SMatrix rotations

* faster, faster

* small update for Euclidean

* bump version and base compat bound

* add one missing test
@mateuszbaran
Copy link
Member

An update: ArrayPartition SA NestedReplacingPowerRepresentation should work in the newest Manifolds.jl. Next up is power group (#500) which addresses this point:

3. How do I do group operations such as, identity_element,vee,hat

@mateuszbaran
Copy link
Member

ArrayPartition SA NestedReplacingPowerRepresentation works up to a small issue with allocation. You have to use X = similar(p) instead of X = allocate(p), and eventually I'm going to let X = allocate(M, p) do the right thing.

@Affie
Copy link
Contributor Author

Affie commented Jul 7, 2022

Thank you very much. I already experimented with it a bit and it seemed to be easy to use. I upgraded to do a solve without errors.

Also, it was a bit faster, thanks!

I made my own identity_element to be able to use StaticArrays with ArrayPartition. What would the recommended way to use different types be. I know about identity_element(M, p) but we don't always have the p available.

I overload a function from DistributedFactorGraphs to start with, something like:

function DFG.getPointIdentity(G::ProductGroup, ::Type{T}=Float64) where T
  M = G.manifold
  return ArrayPartition(map(x->getPointIdentity(x, T), M.manifolds))
end
function DFG.getPointIdentity(G::SpecialOrthogonal{N}, ::Type{T}=Float64) where {N,T}
  return SMatrix{N,N, T}(I)
end

Note: getPointIdentity is originally a way for a user to provide a point type for their own variable type (manifold).

I don't mind adding our defaults with a fallback to identity_element.

@mateuszbaran
Copy link
Member

Also, it was a bit faster, thanks!

Great, that was the idea 🙂 .

I made my own identity_element to be able to use StaticArrays with ArrayPartition. What would the recommended way to use different types be. I know about identity_element(M, p) but we don't always have the p available.

That's a good idea, I will add variants of identity_element that take array eltype as argument. It could similarly work for ArrayPartition vs ProductRepr, I think.

@mateuszbaran
Copy link
Member

To make ArrayPartition work a bit more nicely we also need this: SciML/RecursiveArrayTools.jl#219 .

@mateuszbaran
Copy link
Member

This is an updated test script that uses the new power group and manifold-specific allocation:

using Manifolds
using RecursiveArrayTools: ArrayPartition
using StaticArrays
using Test
using PreallocationTools

G = SpecialEuclidean(2)

# the different options
options = []
#ProductRepr
o = (name = "ProductRepr NestedPowerRepresentation",
     p1 = ProductRepr([0.0,0.0], [1.0 0.0; 0.0 1.0]),
     p2 = ProductRepr([0.0,0.0], [0.0 -1.0; 1.0 0.0]),
     M = PowerManifold(G, NestedPowerRepresentation(), 2));
push!(options, o);

#ArrayPartition
o = (name = "ArrayPartition NestedPowerRepresentation",
     p1 = ArrayPartition([0.0,0.0], [1.0 0.0; 0.0 1.0]),
     p2 = ArrayPartition([0.0,0.0], [0.0 -1.0; 1.0 0.0]),
     M = PowerManifold(G, NestedPowerRepresentation(), 2));
push!(options, o);

#ProductRepr with NestedReplacingPowerRepresentation
o = (name = "ProductRepr NestedReplacingPowerRepresentation",
     p1 = ProductRepr([0.0,0.0], [1.0 0.0; 0.0 1.0]), 
     p2 = ProductRepr([0.0,0.0], [0.0 -1.0; 1.0 0.0]),
     M = PowerManifold(G, NestedReplacingPowerRepresentation(), 2));
push!(options, o);

o = (name = "ProductRepr SA NestedReplacingPowerRepresentation",
     p1 = ProductRepr(SA[0.0,0.0], SA[1.0 0.0; 0.0 1.0]),
     p2 = ProductRepr(SA[0.0,0.0], SA[0.0 -1.0; 1.0 0.0]),
     M = PowerManifold(G, NestedReplacingPowerRepresentation(), 2));
push!(options, o);

#ArrayPartition
o = (name = "ArrayPartition NestedReplacingPowerRepresentation",
     p1 = ArrayPartition([0.0,0.0], [1.0 0.0; 0.0 1.0]),
     p2 = ArrayPartition([0.0,0.0], [0.0 -1.0; 1.0 0.0]),
     M = PowerManifold(G, NestedReplacingPowerRepresentation(), 2));
push!(options, o);

o = (name = "ArrayPartition SA NestedReplacingPowerRepresentation",
     p1 = ArrayPartition(SA[0.0,0.0], SA[1.0 0.0; 0.0 1.0]),
     p2 = ArrayPartition(SA[0.0,0.0], SA[0.0 -1.0; 1.0 0.0]),
     M = PowerManifold(G, NestedReplacingPowerRepresentation(), 2));
push!(options, o);

#ArrayPartition
o = (name = "Vector NestedReplacingPowerRepresentation",
     p1 = [0.0,0.0],
     p2 = [1.0,1.0],
     M = PowerManifold(TranslationGroup(2), NestedReplacingPowerRepresentation(), 2));
push!(options, o);

o = (name = "SVector NestedReplacingPowerRepresentation",
     p1 = SA[0.0,0.0],
     p2 = SA[1.0,1.0],
     M = PowerManifold(TranslationGroup(2), NestedReplacingPowerRepresentation(), 2));
push!(options, o);

## Try to find one that fully Works


for (i,o) in enumerate(options)
  @info "option $i: $(o.name)"
  try
     p = [o.p1,o.p1]
     q = [o.p1,o.p2]
     M = o.M
     log(M, p, q)

     X = allocate(M, p)

     log!(M, X, p, q)
     
     X = [allocate(M, o.p1),allocate(M, o.p1)]
     log!(M, X, p, q)

     exp(M, p, X)

     exp!(M, q, p, X)

     Xc = get_coordinates(M, p, X, DefaultBasis())

     get_coordinates!(M, Xc, p, X, DefaultBasis())

     get_vector(M, p, Xc, DefaultBasis())

     get_vector!(M, X, p, Xc, DefaultBasis())

     MG = PowerGroup(M)

     identity_element(MG)
     Xc2 = vee(MG, p, X)
     hat(MG, p, Xc2)

     # @info "Success of Manifols functions, checking dualcache"
     # dualcache(X)
     @info "Success"
     
     println()
  catch ex
     @error ex
     println()
  end
end

Note the use of allocate(M, p) which is relevant for nested replacing power manifolds, in that is allocates arrays of immutable elements instead of making them mutable as needed for nested non-replacing power manifold. Of course, allocate(M, p) still works for all other cases and is overall the right thing to do.

So the last remaining point here is log(M, Identity(M), p), which just needs to be separately implemented for each group.

@Affie
Copy link
Contributor Author

Affie commented Jul 11, 2022

Thanks, it looks to be working well now and MG = PowerGroup(M) is easy and fits with the rest of the manifolds interface.

@mateuszbaran
Copy link
Member

@Affie I've started working on your point 4 (using Identity instead of identity_element): #503 . What other calls do you need to be able to use Identity in?

@Affie
Copy link
Contributor Author

Affie commented Jul 17, 2022

It looks like log and exp. Perhaps vee and hat also.

@kellertuer
Copy link
Member

For log and exp there are log_lie and exp_lie you want to use, hat and vee are operations on group manifolds, so as far as I can see they work with identity.

@mateuszbaran
Copy link
Member

For log and exp there are log_lie and exp_lie you want to use, hat and vee are operations on group manifolds, so as far as I can see they work with identity.

log_lie(M, q) and log(M, Identity(M), q) are the same only for biinvariant metrics.

@kellertuer
Copy link
Member

kellertuer commented Jul 17, 2022

ah, sure, forgot about that; do you want to use identity(G) and passing on to exp for the others?

@mateuszbaran
Copy link
Member

So I think the current master (soon to be tagged as 0.8.19) has everything listed in this issue implemented 🙂 .

@kellertuer
Copy link
Member

When you register that, keep in mind that some not-so-clever person accidentally increased vision number twice, so the PR for 0.8.19 exists in JuliaRegistry and might be used (again), so it might not automatically merge (I am not sure).

@Affie
Copy link
Contributor Author

Affie commented Jul 22, 2022

Thank you very much!

@Affie Affie closed this as completed Jul 22, 2022
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants