Manifolds.jl
Manifolds.jl copied to clipboard
Respect input array type (static/non-static) unless unreasonable
There are a few places in the codebase where we explicitly created StaticArrays or Arrays regardless of the type of the input.
This isn't great. For one thing, it's always possible the dimension of a manifold is so large that a static array makes no sense. Also, for some of the Christoffel functions, even a 10-dimensional manifold can produce a length 1000 array, which isn't good even if the user did pass a static array.
I propose a few utility function of this form:
@generated maybesize(s::Size{S}) where {S} = prod(S) > 100 ? S : Size{S}()
_similar(x, args...) = similar(x, args...)
_similar(x, s::Size{S}) where {S} = similar(x, S...)
_similar(x::StaticArray, s::Size{S}) where {S} = similar(x, maybesize(s))
_similar(x, ::Type{T}, s::Size{S}) where {S,T} = similar(x, T, S...)
_similar(x::StaticArray, ::Type{T}, s::Size{S}) where {S,T} = similar(x, T, maybesize(s))
_reshape(x, args...) = reshape(x, args...)
_reshape(x, s::Size{S}) where {S} = reshape(x, S...)
_reshape(x::StaticArray, s::Size{S}) where {S} = reshape(x, s)
Unlike the usual similar and reshape, which if Size is passed will always return a sized array, these versions only return a sized array if the original array is sized and if the total length of the requested size is less than some hardcoded cut-off deemed by us to be the maximum size above which it is no longer worth using a sized array.
Whenever size information is known statically (e.g. using the representation_size or manifold_dimension), then one of these functions should be preferred over the original in base.
This makes sense, I'm not fully satisfied with the way we currently handle static arrays.
One thing I'd change in you code is that generated function. Something like this:
@generated maybesize(s::Size{S}) where {S} = prod(S) > 100 ? S : :(s)
looks a bit nicer since the given size object would be reused instead of making a new one.
This change could be coordinated with another change I need to support nested power manifold. When a point is represented by an array of arrays, we should call similar on the innermost array instead of the outermost. The outer array should just be mapped and then the inner arrays would be filled. similar for ProductRepr already does that, it just needs to be extended to arrays of arrays. At the same time we could think about renaming this function as it's not really a standard similar.
PR #98 addresses the _similar part of this issue. The _reshape part needs to be solved separately.