Typed instance-free allocation

NEWS.md

@@ -5,6 +5,12 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [0.15.11] 16/07/2024
+
+### Added
+
+* Function `allocate_as` to generically allocate point and tangent vectors on a manifold without a pre-existing instance but of a particular type.
+
 ## [0.15.10] 19/05/2024
 
 ### Added

Project.toml

@@ -1,7 +1,7 @@
 name = "ManifoldsBase"
 uuid = "3362f125-f0bb-47a3-aa74-596ffd7ef2fb"
 authors = ["Seth Axen <[email protected]>", "Mateusz Baran <[email protected]>", "Ronny Bergmann <[email protected]>", "Antoine Levitt <[email protected]>"]
-version = "0.15.10"
+version = "0.15.11"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

src/ManifoldsBase.jl

@@ -59,7 +59,7 @@ Type `T` is the new number element type [`number_eltype`](@ref), if it is not gi
 the element type of `a` is retained. The `dims` argument can be given for non-nested
 allocation and is forwarded to the function `similar`.
 
-It's behavior can be overriden by a specific manifold, for example power manifold with
+It's behavior can be overridden by a specific manifold, for example power manifold with
 nested replacing representation can decide that `allocate` for `Array{<:SArray}` returns
 another `Array{<:SArray}` instead of `Array{<:MArray}`, as would be done by default.
 """
@@ -86,6 +86,62 @@ function allocate(::AbstractManifold, a, T::Type, dim1::Integer, dims::Integer..
 end
 allocate(::AbstractManifold, a, T::Type, dims::Tuple) = allocate(a, T, dims)
 
+"""
+    allocate_as(M::AbstractManifold, [T:::Type])
+    allocate_as(M::AbstractManifold, F::FiberType, [T:::Type])
+
+Allocate a new point on manifold `M` with optional type given by `T`. Note that `T` is not
+number element type as in [`allocate`](@ref) but rather the type of the entire point to be
+returned.
+
+If `F` is provided, then an element of the corresponding fiber is allocated, assuming it is
+independent of the base point.
+
+To allocate a tangent vector, use ``
+
+# Example
+
+```julia-repl
+julia> using ManifoldsBase
+
+julia> M = ManifoldsBase.DefaultManifold(4)
+DefaultManifold(4; field = ℝ)
+
+julia> allocate_as(M)
+4-element Vector{Float64}:
+ 0.0
+ 0.0
+ 0.0
+ 0.0
+
+julia> allocate_as(M, Array{Float64})
+4-element Vector{Float64}:
+ 0.0
+ 0.0
+ 0.0
+ 0.0
+
+julia> allocate_as(M, TangentSpaceType())
+4-element Vector{Float64}:
+ 0.0
+ 0.0
+ 0.0
+ 0.0
+
+julia> allocate_as(M, TangentSpaceType(), Array{Float64})
+4-element Vector{Float64}:
+ 0.0
+ 0.0
+ 0.0
+ 0.0
+
+```
+"""
+allocate_as(M::AbstractManifold) = similar(Array{Float64}, representation_size(M))
+function allocate_as(M::AbstractManifold, T::Type{<:AbstractArray})
+    return similar(T, representation_size(M))
+end
+
 """
     _pick_basic_allocation_argument(::AbstractManifold, f, x...)
 
@@ -477,7 +533,7 @@ embed!(M::AbstractManifold, Y, p, X) = copyto!(M, Y, p, X)
 
 Embed `p` from manifold `M` an project it back to `M`. For points from `M` this is identity
 but in case embedding is defined for points outside of `M`, this can serve as a way
-to for example remove numerical innacuracies caused by some algorithms.
+to for example remove numerical inaccuracies caused by some algorithms.
 """
 function embed_project(M::AbstractManifold, p)
     return project(M, embed(M, p))
@@ -487,8 +543,8 @@ end
 
 Embed vector `X` tangent at `p` from manifold `M` an project it back to tangent space
 at `p`. For points from that tangent space this is identity but in case embedding is
-defined for tagent vectors from outside of it, this can serve as a way to for example remove
-numerical innacuracies caused by some algorithms.
+defined for tangent vectors from outside of it, this can serve as a way to for example remove
+numerical inaccuracies caused by some algorithms.
 """
 function embed_project(M::AbstractManifold, p, X)
     return project(M, p, embed(M, p, X))
@@ -1219,6 +1275,7 @@ export ×,
     ℝ,
     ℂ,
     allocate,
+    allocate_as,
     angle,
     base_manifold,
     base_point,

src/bases.jl

@@ -280,6 +280,13 @@ const all_uncached_bases{T} = Union{
     DefaultOrthonormalBasis{<:Any,T},
 }
 
+function allocate_as(M::AbstractManifold, ::TangentSpaceType)
+    return similar(Array{Float64}, representation_size(M))
+end
+function allocate_as(M::AbstractManifold, ::TangentSpaceType, T::Type{<:AbstractArray})
+    return similar(T, representation_size(M))
+end
+
 """
     allocate_coordinates(M::AbstractManifold, p, T, n::Int)
 

test/allocation.jl

@@ -67,9 +67,24 @@ ManifoldsBase.representation_size(::AllocManifold3) = (2, 3)
     @test number_eltype(Any[[2.0], [3.0]]) === Float64
     @test number_eltype(typeof([[1.0, 2.0]])) === Float64
 
-    alloc2 = ManifoldsBase.allocate_result(AllocManifold2(), rand)
+    M2 = AllocManifold2()
+    alloc2 = ManifoldsBase.allocate_result(M2, rand)
     @test alloc2 isa Matrix{Float64}
-    @test size(alloc2) == representation_size(AllocManifold2())
+    @test size(alloc2) == representation_size(M2)
 
     @test ManifoldsBase.allocate_result(AllocManifold3(), rand) isa Matrix{ComplexF64}
+
+    an = allocate_as(M2)
+    @test an isa Matrix{Float64}
+    @test size(an) == representation_size(M2)
+    an = allocate_as(M2, Array{Float32})
+    @test an isa Matrix{Float32}
+    @test size(an) == representation_size(M2)
+
+    an = allocate_as(M2, TangentSpaceType())
+    @test an isa Matrix{Float64}
+    @test size(an) == representation_size(M2)
+    an = allocate_as(M2, TangentSpaceType(), Array{Float32})
+    @test an isa Matrix{Float32}
+    @test size(an) == representation_size(M2)
 end