-
Notifications
You must be signed in to change notification settings - Fork 56
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
[WIP] Product manifold #22
Conversation
The canonical product metric is the block diagonal matrix with the individual metrics on the diagonal. I assume you will need to call And how about calling it |
Quite a generic approach. How do you plan to be able to do a |
The main problem with the product metric is at the same time having a default metric structure for products of manifolds with such structure (e.g. I think renaming it to
It's already implemented using reshapes of views ( In FunManifolds, I use an |
src/Product.jl
Outdated
|
||
function det_local_metric(M::MetricManifold{Product{<:Manifold,TRanges,TSizes},ProductMetric}, x) where {TRanges, TSizes} | ||
dets = map(ziptuples(M.manifolds, TRanges, TSizes)) do d | ||
return det_local_metric(d[1], view_element(x, d[2], d[3])) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
If you ziptuples
does what I suppose it's doing, you could take a recursive approach that I have seen in Base
and elsewhere, also to guarantee type inference. Suppose x
was a tuple of coordinate vectors. I guess you're accessing those "tuples" via view_element
.
function det_local_metric(M::MetricManifold{Product{<:Manifold,TRanges,TSizes},ProductMetric}, x) where {TRanges, TSizes}
return det_local_metric(first(M), first(x)) * det_local_metric(tail(M), tail(x))
end
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This won't work right now, x
is not a tuple. view_element
is supposed to take a statically-sized reshaped view of a part of the abstract vector x
.
src/Product.jl
Outdated
end | ||
|
||
function representation_size(M::Product, ::Type{T}) where {T} | ||
return (sum(map(m -> representation_size(m, T), M.manifolds)),) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
As below, this could be
function representation_size(M::Product, ::Type{T}) where {T}
return representation_size(first(M.manifolds), T) + representation_size(tail(M.manifolds), T)
end
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
That would work if M.manifolds
only contains two manifolds, but it could contain an arbitrary number.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Could do a mapreduce
.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
That would work if
M.manifolds
only contains two manifolds, but it could contain an arbitrary number.
Actually the second part recursively upwraps it :) Notice that tail
in the second summand works on all remaining one (again doing the first plus the remaining ones and so on)... nice recursion!
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Oh! Cool, nice approach.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Result of tail
would have to be re-wrapped into a new object of type Product
which doesn't look ideal. mapreduce
makes sense.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Or we would at least need a new function that operates on tuples of manifolds.
Oh, then I have not yet understood your implementation of such views (nor the abbreviation u and su), could you point to a documentation what that does? Otherwise I am looking forward to a documentation of that with an example :) |
test/runtests.jl
Outdated
|
||
for T in types | ||
@testset "Type $T" begin | ||
pts = [convert(T, [1.0, 0.0, 0.0, 0.0, 0.0]), |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
What I meant in my comment is that this construction of a point (just a 5-element array where the first 3 are on S2, the remaining 2 on Rn is quite complicated to read and hard to understand. I can understand that inernally that might be a good representation and that the views (I think I roughly understood them) are a good help, but would prefer a construction (at least additionally) of the type ([1.0, 0.0, 0.0], [0.0, 0.0])
, the elements on the product manifold are ordered anyways (hence the tuple idea).
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I favor this approach too. It also enables us to mix and match representations without hassle and doesn't force us to flatten matrices and such. I also think it's clearer in some cases. For example, SE(n)
's underlying topological manifold is SO(n) \times R(n)
. Its elements can be represented as (n+1) \times (n+1)
matrices which act on elements of R(n+1)
, but since SE(n)
acts on elements of R(n)
, it's cleaner and more efficient to represent its elements as a 2-tuple containing an element of SO(n)
and one of R(n)
.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This could actually be implemented with a custom MPoint
called ProductMPoint
, which stores a list of points and enables dispatch on their types.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
That would be one way, which might not be efficient (though from modelling side I like it). I like how easy @mateuszbaran can compute logs and so on with his internal single-dimension array representation. If a constructor like I proposed and accessorscan be realized (something like p[1]
and p[2]
accessing first and second manifold element of the Product manifold point p
(and that might require a type to dispatch the array access on, maybe). But I would favour that just for ease of use. I am working a lot with product (and even more with power manifolds, even $ \mathcal M^{128\times 128} $ or something, so I really would like to see easy access to single manifolds values.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yeah, I don't know how efficient this approach is. Still working through @mateuszbaran's array-based approach.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'll definitely add a more user-friendly method of creating points on product manifolds once I'm satisfied with the backend.
Tuples are not that great because they are:
- Not
AbstractArray
s, which limits compatibility of the representation with some parts of the Julia ecosystem. - Immutable.
I'd love to have a mutable heterogeneous AbstractArray
but it doesn't exist right now. I'm giving it a serious thought though.
I'll definitely add some convenient accessors.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think you can mutate vectors in a tuple of Vectors/Arrays. Certainly, I'm missing context, but this works.
julia> x = (rand(2), rand(3,3))
([0.178285, 0.857401], [0.525022 0.700751 0.294022; 0.788598 0.253698 0.641965; 0.42019 0.56932 0.875257])
julia> x[1] .= 1
2-element Array{Float64,1}:
1.0
1.0
julia> x
([1.0, 1.0], [0.525022 0.700751 0.294022; 0.788598 0.253698 0.641965; 0.42019 0.56932 0.875257])
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Hmm... it didn't work in FunManifolds because I wanted to be able to represent points as SArray
s but since it's not supported here it might work. I'll try that and see what happens. Things like broadcasting, representation_size
and ambient_space
will be challenging though.
src/Euclidean.jl
Outdated
@@ -19,6 +19,10 @@ struct Euclidean{T<:Tuple} <: Manifold where {T} end | |||
Euclidean(n::Int) = Euclidean{Tuple{n}}() | |||
Euclidean(m::Int, n::Int) = Euclidean{Tuple{m,n}}() | |||
|
|||
function representation_size(::Euclidean{S}, ::Type{T}) where {S,T<:Union{MPoint, TVector, CoTVector}} |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
T
should also permit Array
s here.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Not really, T
is supposed to refer what kind of object is represented (point on the manifold, some element of a vector space like tangent vector, cotangent vector, Riemann tensor etc.). I don't really like this solution, it's not clear what it is, but I thought some functionality of this sort is necessary and this was my first idea.
src/Product.jl
Outdated
end | ||
|
||
function representation_size(M::Product, ::Type{T}) where {T} | ||
return (sum(map(m -> representation_size(m, T), M.manifolds)),) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
That would work if M.manifolds
only contains two manifolds, but it could contain an arbitrary number.
src/Product.jl
Outdated
end | ||
|
||
function representation_size(M::Product, ::Type{T}) where {T} | ||
return (sum(map(m -> representation_size(m, T), M.manifolds)),) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Could do a mapreduce
.
src/Product.jl
Outdated
|
||
function Product(manifolds...) | ||
sizes = map(m -> representation_size(m, MPoint), manifolds) | ||
lengths = map(s -> prod(s), sizes) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Could do map(prod, sizes)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Why wouldn't this implementation of representation_size
work when a product of more than two manifolds is represented?
Thanks for the tips about mapreduce
and prod
. 👍
I was thinking a product of two |
To get the above behavior, we should implement |
I'll document it soon, it's still a WIP and things might change. Especially these two functions. They are supposed to take a statically-sized reshaped views of a part of an array representation of a point, tangent vector etc. |
I've changed the internal representation of product manifolds quite a lot. Finally, it's a mix of a few approached discussed here: points, vectors etc. are represented by objects of type What do you think about the new approach? |
I think the new approach is really cool (thanks for setting up the benchmark suite, by the way. Awesome!), and I'm kind of shocked that it performs so well. I'm still trying to understand how it works. How would it handle a mixture of vector and matrix manifolds? I'm not certain what to pass to Currently, julia> R2 = Manifolds.Euclidean(2);
julia> M = Manifolds.ProductManifold(R2, R2);
julia> x = randn(4)
4-element Array{Float64,1}:
0.7683675344912058
0.045762904704441366
1.281251818793342
1.148563618371174
julia> P = Manifolds.ProductView(M, x)
4-element Manifolds.ProductView{Manifolds.ProductManifold{Tuple{Euclidean{Tuple{2}},Euclidean{Tuple{2}}},(1:2, 3:4),Tuple{Tuple{2},Tuple{2}}},Float64,1,Array{Float64,1},Tuple{Manifolds.SizedAbstractArray{Tuple{2},Float64,1,1,SubArray{Float64,1,Array{Float64,1},Tuple{UnitRange{Int64}},true}},Manifolds.SizedAbstractArray{Tuple{2},Float64,1,1,SubArray{Float64,1,Array{Float64,1},Tuple{UnitRange{Int64}},true}}}}:
0.7683675344912058
0.045762904704441366
1.281251818793342
1.148563618371174
julia> M = Manifolds.ProductManifold(R2, R2, R2);
julia> x = randn(6)
6-element Array{Float64,1}:
0.6630836648853156
0.2331471356770406
0.2964616068062833
-0.14207468080177554
-1.5358504686848118
0.3998359958877595
julia> P = Manifolds.ProductView(M, x)
ERROR: MethodError: no method matching ziptuples(::Tuple{UnitRange{Int64},UnitRange{Int64},UnitRange{Int64}}, ::Type{Tuple{Tuple{2},Tuple{2},Tuple{2}}})
Closest candidates are:
ziptuples(::Tuple{Vararg{Any,N}}, ::Tuple{Vararg{Any,M}}) where {N, M} at /Users/saxen/projects/Manifolds/src/utils.jl:9
ziptuples(::Tuple{Vararg{Any,N}}, ::Tuple{Vararg{Any,M}}, ::Tuple{Vararg{Any,L}}) where {N, M, L} at /Users/saxen/projects/Manifolds/src/utils.jl:23
ziptuples(::Tuple{Vararg{Any,N}}, ::Tuple{Vararg{Any,M}}, ::Tuple{Vararg{Any,L}}, ::Tuple{Vararg{Any,K}}) where {N, M, L, K} at /Users/saxen/projects/Manifolds/src/utils.jl:37
Stacktrace:
[1] Manifolds.ProductView(::Manifolds.ProductManifold{Tuple{Euclidean{Tuple{2}},Euclidean{Tuple{2}},Euclidean{Tuple{2}}},(1:2, 3:4, 5:6),Tuple{Tuple{2},Tuple{2},Tuple{2}}}, ::Array{Float64,1}) at /Users/saxen/projects/Manifolds/src/ProductManifold.jl:50
[2] top-level scope at none:0 |
I'm glad you like the new approach 😃. I don't really know why it's so fast either, I was copying pieces of Julia's I've added a product of three manifolds including a matrix one (SO(2)) to tests. I've added I didn't really think yet about how |
It caused performance regressions in some run for ProductManifold with ProductArray. I have no idea why, but removing tune.json and re-running seems to help. At least it has helped once.
I've restructured things a bit to allow products of arbitrary manifolds (not necessarily array-based). It has surprisingly good performance, generally on par with |
Codecov Report
@@ Coverage Diff @@
## master #22 +/- ##
=========================================
Coverage ? 66.84%
=========================================
Files ? 10
Lines ? 555
Branches ? 0
=========================================
Hits ? 371
Misses ? 184
Partials ? 0
Continue to review full report at Codecov.
|
Codecov Report
@@ Coverage Diff @@
## master #22 +/- ##
=========================================
Coverage ? 72.58%
=========================================
Files ? 10
Lines ? 591
Branches ? 0
=========================================
Hits ? 429
Misses ? 162
Partials ? 0
Continue to review full report at Codecov.
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
That's in total quite some progress! Thanks. Looking forward to powermanifold and tangent bundle in this efficient form!
""" | ||
function proj_product(x::ProductArray, i::Integer) | ||
return copy(x.parts[i]) | ||
end |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I find the function name a little strange – but I think you discussed that already? For me its not a project but merely an access (but with copy). Btw. Why do you copy? Without that one could just overload [I]
access for ProductArray
?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think projections were discussed in the context of embedded manifolds. I've seen such functions from the Cartesian product to one of its factors also called projections but we can use a different name here. Unfortunately, getindex
/setindex!
(or [I]
) is not an option, they are already defined in a different way, as required by broadcasting. I could probably write much more complex broadcasting rules that don't require getindex
and setindex!
defined this way but I don't think it's worth the trouble.
I'll remove that copy
, I was experimenting with some things and didn't remove it.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Ah I did not see that difference and was just thinking about ease of use, since (just as a user view again) having a point p
on a Productmanifold
a just natural thing to access the i
th entry would be p[I]
; but yes I am also not yet so familiar with broadcast rules.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
So, how do we name this function then?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Well, it's a submanifold of the original PowerManifold
, so maybe submanifold
could return a specific manifold from a ProductManifold
and submanifold_point
could return the point? The downside to this name is that one would expect to be able to provide any number of indices to build a submanifold, not just a single one.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I prefer if we're going to use submanifold
to get a subset of manifolds that we use something like submanifold_
here. I think submanifold_point
could still work, because no manifold is provided, and a TVector
and CoTVector
are points on manifolds in their own rights, but since we haven't been using "point" in this general way so far, I can see how it would be confusing. I think _component
works. If this and power manifold is the only place where we'll be using these functions, it might make sense to use factor_
instead of submanifold_
, since that is the more descriptive term.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I don't know, I'm not that good at naming things. I think I'll just let you and @kellertuer pick names for this and prod_point
.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Okay, here're my thoughts. Let's keep submanifold
, since we may have more uses for that function later. I don't think we should name this submanifold_point
, since one would expect that to return something similar to the original ProductArray
. Instead, it should be submanifold_component(s)
, since the user gets a container of parts. I don't know whether the function name should be plural.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
OK, I've renamed it to submanifold_component
. Singular seems more appropriate, right now components can be accessed only one at a time since there is no obvious way to handle multi-component access for ProductArray
(new ProductArray
? a tuple with views? an array with views?)
BTW, how do we rename product_point
now?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yay! Picking names! Supercalifragilistic_expialigetic
(sorry its late). I was liking the name access(p,index)
or access(p,indices)
(which would also get in handy if we figure out how to rename that to getindex and get really p[i]
-things. Otherwise I am also fine with manifold_component
without plural, since for a single index its a point on a manifold and even for a container of parts, that could be seen as a point on a submanifold that is still a product manifold; the consistent way for both is hence the singular.
|
||
function proj_product(x::ProductTVector, i::Integer) | ||
return x.parts[i] | ||
end |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
similar as above – and why do you then not copy here?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Let's discuss it above and I'll change both methods to make them consistent.
tested. | ||
- `retraction_methods = []`: retraction methods that will be tested. | ||
""" | ||
function test_manifold(M::Manifold, pts::AbstractVector; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
It's great to have this in such generality, I think manifold_dimension
and injectivity_radius
are missing, but might also be done one by one per manifold. Similarly is_manifold_point
(also testing the errors that might occur) is missing and is_tangent_vector
.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
and great that we have codecov running 👍
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'd like to have the whole interface for manifolds tested here and so I expand it slowly. manifold_dimension
and injectivity_radius
are there, is_manifold_point
and is_tangent_vector
are partially tested (only the positive case). I think the negative case (when an error is thrown) should be tested separately for each manifold.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yes, I also think, that the negative cases can only be tested individually by carefully constructing invalid cases (if they even exist) and check for the right error messages.
Codecov Report
@@ Coverage Diff @@
## master #22 +/- ##
=========================================
Coverage ? 72.95%
=========================================
Files ? 10
Lines ? 599
Branches ? 0
=========================================
Hits ? 437
Misses ? 162
Partials ? 0
Continue to review full report at Codecov.
|
Next, the desired point on the product manifold can be obtained by calling | ||
`Manifolds.prod_point(Mshape, [1.0, 0.0, 0.0], [-3.0, 2.0])`. | ||
""" | ||
function prod_point(M::ShapeSpecification, pts...) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Would be nice to bring this name in line with submanifold_component
or whatever we decide, since this isn't unique to points.
Building a julia> E = Euclidean(3)
Euclidean{Tuple{3}}()
julia> P = Manifolds.ProductManifold(E,E,E,E)
Manifolds.ProductManifold{NTuple{4,Euclidean{Tuple{3}}}}((Euclidean{Tuple{3}}(), Euclidean{Tuple{3}}(), Euclidean{Tuple{3}}(), Euclidean{Tuple{3}}()))
julia> shape = Manifolds.ShapeSpecification(P.manifolds...)
Manifolds.ShapeSpecification{(1:3, 4:6, 7:9, 10:12),NTuple{4,Tuple{3}}}()
julia> x = randn(12);
julia> a = Manifolds.ProductArray(shape, x)
ERROR: MethodError: no method matching Manifolds.ProductArray(::Type{Manifolds.ShapeSpecification{(1:3, 4:6, 7:9, 10:12),NTuple{4,Tuple{3}}}}, ::Array{Float64,1})
Closest candidates are:
Manifolds.ProductArray(::Type{Manifolds.ShapeSpecification{TRanges,Tuple{Size1,Size2}}}, ::TData<:AbstractArray{T,N}) where {TRanges, Size1, Size2, T, N, TData<:AbstractArray{T,N}} at /Users/saxen/projects/Manifolds/src/ProductManifold.jl:118
Manifolds.ProductArray(::Type{Manifolds.ShapeSpecification{TRanges,Tuple{Size1,Size2,Size3}}}, ::TData<:AbstractArray{T,N}) where {TRanges, Size1, Size2, Size3, T, N, TData<:AbstractArray{T,N}} at /Users/saxen/projects/Manifolds/src/ProductManifold.jl:124
Manifolds.ProductArray(::Manifolds.ShapeSpecification{TRanges,TSizes}, ::TData<:AbstractArray{T,N}) where {TM, TRanges, TSizes, T, N, TData<:AbstractArray{T,N}} at /Users/saxen/projects/Manifolds/src/ProductManifold.jl:131
Stacktrace:
[1] Manifolds.ProductArray(::Manifolds.ShapeSpecification{(1:3, 4:6, 7:9, 10:12),NTuple{4,Tuple{3}}}, ::Array{Float64,1}) at /Users/saxen/projects/Manifolds/src/ProductManifold.jl:131
[2] top-level scope at none:0``` |
OK, I've fixed that bug with products of more than three manifolds, re-exported |
Will |
I've exported |
I've cleaned up the interface a bit. I'd like to merge it and resolve remaining issues in separate PRs, it'll only get harder and harder to merge other PRs with the changes made here. |
I've started working on the product manifold. Functionality is going quite nicely, remaining challenges:
zero_tangent_vector!
for sphere is an interesting example.