Nested power representation with replacement

[mateuszbaran]

That's the thing I've been recently talking about, i.e. a nested power manifold that doesn't require mutable representation of the wrapped manifold.

[codecov]

Codecov Report

Merging #73 (18d07a9) into master (194e27e) will increase coverage by 0.15%.
The diff coverage is 100.00%.

❗ Current head 18d07a9 differs from pull request most recent head 5bfecd6. Consider uploading reports for the commit 5bfecd6 to get more accurate results

@@            Coverage Diff             @@
##           master      #73      +/-   ##
==========================================
+ Coverage   99.70%   99.85%   +0.15%     
==========================================
  Files          14       14              
  Lines        1351     1426      +75     
==========================================
+ Hits         1347     1424      +77     
+ Misses          4        2       -2     

Impacted Files	Coverage Δ
src/ManifoldsBase.jl	99.08% <ø> (ø)
src/DefaultManifold.jl	100.00% <100.00%> (ø)
src/PowerManifold.jl	100.00% <100.00%> (+0.85%)	⬆️
src/bases.jl	100.00% <100.00%> (ø)

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 194e27e...5bfecd6. Read the comment docs.

[kellertuer]

This looks interesting, though the name might be irritating due to its length; maybe it would help to have a concrete example how these two types of nested actions differ?

[mateuszbaran]

The goal is to have a nested representation that supports immutable representation of points and tangent vectors on the wrapped manifold, for example SE(2)^n where a point on this power manifold is stored in-line in a vector. I'll make a more comprehensive example later.

Also, I'm definitely open to changing the name but so far I haven't come up with anything better.

[kellertuer]

No I think the name sounds fine, if (but that would be breaking) I would call the other now NestedMutatingPowerRepresentation, maybe even with a common abstract super type.

Just for users to start with giving the details on similarities and differences would be neat.

[mateuszbaran]

I'm not quite sure how to phrase it. Here is an example of a power manifold of SE(2) with in-line storage (all point on SE(2) are stored in a continuous block of memory):

julia> using Manifolds, ManifoldsBase, StaticArrays

julia> R2 = Rotations(2)
Rotations(2)

julia> G = SpecialEuclidean(2)
SpecialEuclidean(2)

julia> N = 5

julia> GN = PowerManifold(G, ManifoldsBase.NestedReplacingPowerRepresentation(), N)
PowerManifold(SpecialEuclidean(2), ManifoldsBase.NestedReplacingPowerRepresentation(), N)

julia> q = [1.0 0.0; 0.0 1.0]

julia> p1 = [ProductRepr(SVector{2,Float64}([i - 0.1, -i]), SMatrix{2,2,Float64}(exp(R2, q, hat(R2, q, i)))) for i in 1:N]
5-element Vector{ProductRepr{Tuple{SVector{2, Float64}, SMatrix{2, 2, Float64, 4}}}}:
 ProductRepr{Tuple{SVector{2, Float64}, SMatrix{2, 2, Float64, 4}}}(([0.9, -1.0], [0.5403023058681398 -0.8414709848078965; 0.8414709848078965 0.5403023058681398]))
 ProductRepr{Tuple{SVector{2, Float64}, SMatrix{2, 2, Float64, 4}}}(([1.9, -2.0], [-0.4161468365471424 -0.9092974268256817; 0.9092974268256817 -0.4161468365471424]))
 ProductRepr{Tuple{SVector{2, Float64}, SMatrix{2, 2, Float64, 4}}}(([2.9, -3.0], [-0.9899924966004454 -0.1411200080598672; 0.1411200080598672 -0.9899924966004454]))
 ProductRepr{Tuple{SVector{2, Float64}, SMatrix{2, 2, Float64, 4}}}(([3.9, -4.0], [-0.6536436208636119 0.7568024953079282; -0.7568024953079282 -0.6536436208636119]))
 ProductRepr{Tuple{SVector{2, Float64}, SMatrix{2, 2, Float64, 4}}}(([4.9, -5.0], [0.28366218546322625 0.9589242746631385; -0.9589242746631385 0.28366218546322625]))

julia> p2 = [ProductRepr(SVector{2,Float64}([i - 0.1, -i]), SMatrix{2,2,Float64}(exp(R2, q, hat(R2, q, -i)))) for i in 1:N]
5-element Vector{ProductRepr{Tuple{SVector{2, Float64}, SMatrix{2, 2, Float64, 4}}}}:
 ProductRepr{Tuple{SVector{2, Float64}, SMatrix{2, 2, Float64, 4}}}(([0.9, -1.0], [0.5403023058681398 0.8414709848078965; -0.8414709848078965 0.5403023058681398]))
 ProductRepr{Tuple{SVector{2, Float64}, SMatrix{2, 2, Float64, 4}}}(([1.9, -2.0], [-0.4161468365471424 0.9092974268256817; -0.9092974268256817 -0.4161468365471424]))
 ProductRepr{Tuple{SVector{2, Float64}, SMatrix{2, 2, Float64, 4}}}(([2.9, -3.0], [-0.9899924966004454 0.1411200080598672; -0.1411200080598672 -0.9899924966004454]))
 ProductRepr{Tuple{SVector{2, Float64}, SMatrix{2, 2, Float64, 4}}}(([3.9, -4.0], [-0.6536436208636119 -0.7568024953079282; 0.7568024953079282 -0.6536436208636119]))
 ProductRepr{Tuple{SVector{2, Float64}, SMatrix{2, 2, Float64, 4}}}(([4.9, -5.0], [0.28366218546322625 -0.9589242746631385; 0.9589242746631385 0.28366218546322625]))

julia> X = similar(p1)
5-element Vector{ProductRepr{Tuple{SVector{2, Float64}, SMatrix{2, 2, Float64, 4}}}}:
 ProductRepr{Tuple{SVector{2, Float64}, SMatrix{2, 2, Float64, 4}}}(([0.0, 0.0], [0.0 0.0; 0.0 0.0]))
 ProductRepr{Tuple{SVector{2, Float64}, SMatrix{2, 2, Float64, 4}}}(([0.0, 0.0], [0.0 0.0; 0.0 0.0]))
 ProductRepr{Tuple{SVector{2, Float64}, SMatrix{2, 2, Float64, 4}}}(([0.0, 0.0], [0.0 0.0; 0.0 0.0]))
 ProductRepr{Tuple{SVector{2, Float64}, SMatrix{2, 2, Float64, 4}}}(([0.0, 0.0], [0.0 0.0; 0.0 0.0]))
 ProductRepr{Tuple{SVector{2, Float64}, SMatrix{2, 2, Float64, 4}}}(([0.0, 0.0], [0.0 0.0; 0.0 0.0]))

julia> log!(GN, X, p1, p2)
5-element Vector{ProductRepr{Tuple{SVector{2, Float64}, SMatrix{2, 2, Float64, 4}}}}:
 ProductRepr{Tuple{SVector{2, Float64}, SMatrix{2, 2, Float64, 4}}}(([0.0, 0.0], [0.0 2.0; -2.0 0.0]))
 ProductRepr{Tuple{SVector{2, Float64}, SMatrix{2, 2, Float64, 4}}}(([0.0, 0.0], [0.0 -2.2831853071795862; 2.2831853071795862 0.0]))
 ProductRepr{Tuple{SVector{2, Float64}, SMatrix{2, 2, Float64, 4}}}(([0.0, 0.0], [0.0 -0.28318530717958645; 0.28318530717958645 0.0]))
 ProductRepr{Tuple{SVector{2, Float64}, SMatrix{2, 2, Float64, 4}}}(([0.0, 0.0], [0.0 1.7168146928204135; -1.7168146928204135 0.0]))
 ProductRepr{Tuple{SVector{2, Float64}, SMatrix{2, 2, Float64, 4}}}(([0.0, 0.0], [0.0 -2.566370614359173; 2.566370614359173 0.0]))

I guess there should be perhaps a centralized place that explains strengths and weaknesses of each representation? Maybe in https://juliamanifolds.github.io/Manifolds.jl/stable/manifolds/power.html ?

[kellertuer]

Thanks for the example, I think I have an idea of what the difference is. Yes that should be documented (all 3 types atcually) at https://juliamanifolds.github.io/Manifolds.jl/stable/manifolds/power.html

Let me try to phrase it in my own words and we'll see whether this is correct.
Let $p\in \mathcal M$. be a point, for simplicity just an array of length $m$ and $P\in \mathcal M^n$.

1. ArrayRepresentation: This is kind of the easy or maybe also called flat-array-approach: The data on the power manifold is stored as one array. Here $P$ is a matrix of m rows, n columns. This array can be a classic array, a hybrid or a static one. Here, mutating variants work for the classical cases and the non.mutating (hybrid/static) cases allocate new memory.

2. NestedPowerRepresentation: We store $P$ as a vector of length $n$ containing vectors, i.e. P = [p1, p2, p3,..., pn]. For classical arrays this would mutate the nested arrays. The outer one should always stay the same anywaysand just the inner oner are mutated, right? For nested of nested ones (power of power for example) even the most inner are mutated only, right?

3. The new case NestedReplacingPowerRepresentation, we would store $P$ also as a vector of elements, but these might be static arrays or ProductRepr, so non-mutating ones. Then we force a function here to not mutate on the inner level but replace. This is similar to using 1. with a static variant where also always the whole data is replaced.

In the beginning (this morning, and when starting to write this), I thought whether we would need a depth level for 3, but we don't, since for deeper nested structures one would have a PowerManifold(PowerManifold(...)) anyways and then each Powermanifold can specify their Nested approach. Let's look at a few cases to see whether I understood.

• A nested power (outer) of an array power (inner) where the inner one is static has to be a Replacing one
• an array power (outer) of a nested power (inner) can not happen?
• a nested power (outer,A) of a nested inner (inner,B) we can have
  • A, B as in 2, if we use classic arrays (also for the manifold points inner
  • A,B as in 3, if the manifold has a type that is non-mutating
  • A as in 2, B as in 3 – not sure if this would work for non-mutating manifold points? If not this might just error at the first try of mutation and I'd be fine with that.
  • A as in 3, B as in 2 – for mutating points (arras&Sphere) on the manifold – would just be as the first case with performance reduction due to more copying in-between?

And maybe this nested of nested stuff is already too technical to be mentioned in the docs.

Edit: For your last 2 lines I am not 100% sure, but allocating X when it gets replaced afterwards anyways is not that nice performance wise, but makes these types (non mutating inner ones) work with our general framework, right?

[mateuszbaran]

1. ArrayRepresentation: This is kind of the easy or maybe also called flat-array-approach: The data on the power manifold is stored as one array. Here $P$ is a matrix of m rows, n columns. This array can be a classic array, a hybrid or a static one. Here, mutating variants work for the classical cases and the non.mutating (hybrid/static) cases allocate new memory.

Yes, that's right.

2. NestedPowerRepresentation: We store $P$ as a vector of length $n$ containing vectors, i.e. P = [p1, p2, p3,..., pn]. For classical arrays this would mutate the nested arrays. The outer one should always stay the same anywaysand just the inner oner are mutated, right? For nested of nested ones (power of power for example) even the most inner are mutated only, right?

That's also right, except IIRC we currently only support power-of-power by adding more indices, so no true nesting.

3. The new case NestedReplacingPowerRepresentation, we would store $P$ also as a vector of elements, but these might be static arrays or ProductRepr, so non-mutating ones. Then we force a function here to not mutate on the inner level but replace. This is similar to using 1. with a static variant where also always the whole data is replaced.

Yes, exactly.

• A nested power (outer) of an array power (inner) where the inner one is static has to be a Replacing one

Do we... actually want to support that?

• an array power (outer) of a nested power (inner) can not happen?

Yes, outside of degenerate cases like a zero-index power manifold, it can't happen.

a nested power (outer,A) of a nested inner (inner,B) we can have

• A, B as in 2, if we use classic arrays (also for the manifold points inner
• A,B as in 3, if the manifold has a type that is non-mutating
• A as in 2, B as in 3 – not sure if this would work for non-mutating manifold points? If not this might just error at the first try of mutation and I'd be fine with that.
• A as in 3, B as in 2 – for mutating points (arras&Sphere) on the manifold – would just be as the first case with performance reduction due to more copying in-between?

I'll have to think about it a bit more, it's not something I expected to support.

Edit: For your last 2 lines I am not 100% sure, but allocating X when it gets replaced afterwards anyways is not that nice performance wise, but makes these types (non mutating inner ones) work with our general framework, right?

It's just an example. Sometimes you need to overwrite X multiple times. You surely can use non-mutating log but the point I'm trying to make here is that storage of X is efficient and, with a few tweaks, log!(GN, X, p1, p2) would be non-allocating.

[kellertuer]

That's also right, except IIRC we currently only support power-of-power by adding more indices, so no true nesting.

If you follow the (M^n)^m operator definitions, that is the same as M^(n,m) yes, so the easy generators avoid nestedness. So for the usual user, this is the case, even when you use the constructor PowerManifold(PowerManifold(M,3),3). Then the outer one inherits the representation type of the inner. But you can trick this (and I think that is ok) if you do PowerManifold(PowerManifold(M, ArrayPowerRepresentation(),3), NestedRepresentation(),4). So specifying the outer nested type does true nesting (even if you specify the outer the same to the inner, since specifying it calls the generic AbstractManifold case.

3. [...]

Yes, exactly.

Cool. Then I understood this new case!

• A nested power (outer) of an array power (inner) where the inner one is static has to be a Replacing one

Do we... actually want to support that?

It is possible at least (see above), whether it is useful or we recommend using that is a different thing.

a nested power (outer,A) of a nested inner (inner,B) we can have

• A, B as in 2, if we use classic arrays (also for the manifold points inner
• A,B as in 3, if the manifold has a type that is non-mutating
• A as in 2, B as in 3 – not sure if this would work for non-mutating manifold points? If not this might just error at the first try of mutation and I'd be fine with that.
• A as in 3, B as in 2 – for mutating points (arras&Sphere) on the manifold – would just be as the first case with performance reduction due to more copying in-between?

I'll have to think about it a bit more, it's not something I expected to support.

The first two are default ones, which should work (though see atop one has to trick a little to get them actually nested and not just stacked). For the third case – we should not support that, for the fourth – it should work but it'll be slow, and there is no real use in this I think.

Edit: For your last 2 lines I am not 100% sure, but allocating X when it gets replaced afterwards anyways is not that nice performance wise, but makes these types (non mutating inner ones) work with our general framework, right?

It's just an example. Sometimes you need to overwrite X multiple times. You surely can use non-mutating log but the point I'm trying to make here is that storage of X is efficient and, with a few tweaks, log!(GN, X, p1, p2) would be non-allocating.

Cool. My main point would be, that most generic things I write would use log! (and not log) and with this new representation type, the static things you use for the inner manifold would still work, that is neat! If they even would be non-allocating – even better!

[mateuszbaran]

If you follow the (M^n)^m operator definitions, that is the same as M^(n,m) yes, so the easy generators avoid nestedness. So for the usual user, this is the case, even when you use the constructor PowerManifold(PowerManifold(M,3),3). Then the outer one inherits the representation type of the inner. But you can trick this (and I think that is ok) if you do PowerManifold(PowerManifold(M, ArrayPowerRepresentation(),3), NestedRepresentation(),4). So specifying the outer nested type does true nesting (even if you specify the outer the same to the inner, since specifying it calls the generic AbstractManifold case.

Right, I just didn't think of that.

The first two are default ones, which should work (though see atop one has to trick a little to get them actually nested and not just stacked). For the third case – we should not support that, for the fourth – it should work but it'll be slow, and there is no real use in this I think.

I've thought a bit about it, and maybe I'm missing something, but I've summarized my thoughts in the table below.

Outer (rows) \ Inner (columns)	Array	Nested modifying	Nested replacing
Array	It should work but but I don't see any benefits over single-level array	it may work in some cases but when it does, it effectively works like nested modifying with more indices	it may work in some cases but when it does, it effectively works like nested replacing with more indices
Nested modifying	It should work but why not just using Array power manifold for all indices?	It should work but but I don't see any benefits over single-level nested modifying	this does actually make some sense, f.e. if someone needs to work with things like Vector{Vector{<:ProductRepr}} This variant should correctly handle it; still, Matrix{<:ProductRepr} with nested replacing would be preferable unless someone absolutely needs to use nested vectors.
Nested replacing	It should work but why not just using Array power manifold for all indices?	yes, it's effectively a slower double-level nested modifying	It should work but but I don't see any benefits over single-level nested replacing

In essence, I'd just recommend considering adding more indices to the inner power manifold. Other combinations should work, at least in some cases. I wouldn't add tests for them to our test suite unless some user specifically asks for it.

Cool. My main point would be, that most generic things I write would use log! (and not log) and with this new representation type, the static things you use for the inner manifold would still work, that is neat! If they even would be non-allocating – even better!

That's right, modifying variants still take priority over the non-modifying ones but I'm actually open to adding some non-mutating ones alongside mutating ones, if someone needs that. I expect Euclidean, Rotations and ProductManifold will need separate non-mutating implementations to offer better performance for power manifold of SE(n), since it looks like it would likely be desirable f.e. for that: JuliaRobotics/DistributedFactorGraphs.jl#763 .

[kellertuer]

Thanks for the table, yes, that's a good overview and I agree, that None of these are that good, and the default (adding dimenions to one PowerManifold) should stay the default. If one chooses one of these and goes the extra mile of construction these – the user will have a good reason.

[kellertuer]

Could you maybe check whether you can provide a nice test for the allocation of the result of get_coordinates for both an abstract basis (maybe DefaultOrthogonalBasis and a CachedBasis?

[mateuszbaran]

The AbstractBasis case is covered, and the CachedBasis case doesn't actually seem needed.

[mateuszbaran]

That is, it doesn't seem to be needed for anything other than not having a harmless ambiguity.

[kellertuer]

You're right, it's twice the CachedBasis, I misread the Codecov result. I am not sure whether this is harmless, but I feel if we have the code, we should have a test?

[mateuszbaran]

OK, I've added a test.

[mateuszbaran]

I'm not sure if I should export the new representation? And are you OK with the current naming?

[kellertuer]

I think we should, since everybody should be able to use it.

I think the name is a little long, but I do not have a better name in mind at the moment.