Typed instance-free allocation

.github/workflows/ci.yml

@@ -11,7 +11,7 @@ jobs:
     runs-on: ${{ matrix.os }}
     strategy:
       matrix:
-        julia-version: ["1.6", "1.9", "~1.10.0-0"]
+        julia-version: ["1.6", "1.10", "~1.11.0-0"]
         os: [ubuntu-latest, macOS-latest, windows-latest]
     steps:
       - uses: actions/checkout@v4

.github/workflows/documenter.yml

@@ -13,10 +13,10 @@ jobs:
       - uses: actions/checkout@v4
       - uses: quarto-dev/quarto-actions/setup@v2
         with:
-          version: 1.3.353
+          version: "1.5.55"
       - uses: julia-actions/setup-julia@latest
         with:
-          version: 1.9
+          version: "1.10"
       - name: Julia Cache
         uses: julia-actions/cache@v2
       - name: Cache Quarto

.github/workflows/format.yml

@@ -13,7 +13,7 @@ jobs:
       - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v2
         with:
-          version: 1.9
+          version: "1.10"
       - name: Install JuliaFormatter and format
         run: |
           julia  -e 'using Pkg; Pkg.add(PackageSpec(name="JuliaFormatter"))'

NEWS.md

@@ -5,6 +5,18 @@ 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] 28/07/2024
+
+### Added
+
+* Function `allocate_on` to generically allocate point and tangent vectors on a manifold without a pre-existing instance but of a particular type.
+* Function `default_type` to get the default type of points and tangent vectors for a manifold.
+* Package extension for the `Quaternions.jl` package that handles allocation.
+
+### Changed
+
+* Default allocation method was made more robust to custom promotion functions.
+
 ## [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"
@@ -12,11 +12,13 @@ Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 
 [weakdeps]
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+Quaternions = "94ee1d12-ae83-5a48-8b1c-48b8ff168ae0"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [extensions]
 ManifoldsBasePlotsExt = "Plots"
+ManifoldsBaseQuaternionsExt = "Quaternions"
 ManifoldsBaseRecursiveArrayToolsExt = "RecursiveArrayTools"
 ManifoldsBaseStatisticsExt = "Statistics"
 
@@ -35,6 +37,7 @@ Requires = "1"
 ReverseDiff = "1"
 StaticArrays = "1"
 Statistics = "1.6"
+Quaternions = "0.7"
 Test = "1.6"
 julia = "1.6"
 
@@ -44,11 +47,12 @@ DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+Quaternions = "94ee1d12-ae83-5a48-8b1c-48b8ff168ae0"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "Aqua", "DoubleFloats", "ForwardDiff", "OrdinaryDiffEq", "Plots", "ReverseDiff", "StaticArrays", "Statistics", "RecursiveArrayTools"]
+test = ["Test", "Aqua", "DoubleFloats", "ForwardDiff", "OrdinaryDiffEq", "Plots", "Quaternions", "ReverseDiff", "StaticArrays", "Statistics", "RecursiveArrayTools"]

ext/ManifoldsBaseQuaternionsExt.jl

@@ -0,0 +1,23 @@
+module ManifoldsBaseQuaternionsExt
+
+if isdefined(Base, :get_extension)
+    using ManifoldsBase
+    using ManifoldsBase: ℍ
+    using Quaternions
+else
+    # imports need to be relative for Requires.jl-based workflows:
+    # https://github.com/JuliaArrays/ArrayInterface.jl/pull/387
+    using ..ManifoldsBase
+    using ..ManifoldsBase: ℍ
+    using ..Quaternions
+end
+
+@inline function ManifoldsBase.allocate_result_type(
+    ::AbstractManifold{ℍ},
+    f::TF,
+    args::Tuple{},
+) where {TF}
+    return QuaternionF64
+end
+
+end

ext/ManifoldsBaseRecursiveArrayToolsExt/ManifoldsBaseRecursiveArrayToolsExt.jl

@@ -9,6 +9,7 @@ if isdefined(Base, :get_extension)
 
     import ManifoldsBase:
         allocate,
+        allocate_on,
         allocate_result,
         default_inverse_retraction_method,
         default_retraction_method,
@@ -40,6 +41,7 @@ else
 
     import ..ManifoldsBase:
         allocate,
+        allocate_on,
         allocate_result,
         default_inverse_retraction_method,
         default_retraction_method,

ext/ManifoldsBaseRecursiveArrayToolsExt/ProductManifoldRecursiveArrayToolsExt.jl

@@ -5,6 +5,23 @@ function allocate(x::ArrayPartition, T::Type)
     return ArrayPartition(map(t -> allocate(t, T), submanifold_components(x))...)
 end
 
+allocate_on(M::ProductManifold) = ArrayPartition(map(N -> allocate_on(N), M.manifolds)...)
+function allocate_on(M::ProductManifold, ::Type{ArrayPartition{T,U}}) where {T,U}
+    return ArrayPartition(map((N, V) -> allocate_on(N, V), M.manifolds, U.parameters)...)
+end
+function allocate_on(M::ProductManifold, ft::TangentSpaceType)
+    return ArrayPartition(map(N -> allocate_on(N, ft), M.manifolds)...)
+end
+function allocate_on(
+    M::ProductManifold,
+    ft::TangentSpaceType,
+    ::Type{ArrayPartition{T,U}},
+) where {T,U}
+    return ArrayPartition(
+        map((N, V) -> allocate_on(N, ft, V), M.manifolds, U.parameters)...,
+    )
+end
+
 @inline function allocate_result(M::ProductManifold, f)
     return ArrayPartition(map(N -> allocate_result(N, f), M.manifolds))
 end

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_on(M::AbstractManifold, [T:::Type])
+    allocate_on(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_on(M)
+4-element Vector{Float64}:
+ 0.0
+ 0.0
+ 0.0
+ 0.0
+
+julia> allocate_on(M, Array{Float64})
+4-element Vector{Float64}:
+ 0.0
+ 0.0
+ 0.0
+ 0.0
+
+julia> allocate_on(M, TangentSpaceType())
+4-element Vector{Float64}:
+ 0.0
+ 0.0
+ 0.0
+ 0.0
+
+julia> allocate_on(M, TangentSpaceType(), Array{Float64})
+4-element Vector{Float64}:
+ 0.0
+ 0.0
+ 0.0
+ 0.0
+
+```
+"""
+allocate_on(M::AbstractManifold) = similar(Array{Float64}, representation_size(M))
+function allocate_on(M::AbstractManifold, T::Type{<:AbstractArray})
+    return similar(T, representation_size(M))
+end
+
 """
     _pick_basic_allocation_argument(::AbstractManifold, f, x...)
 
@@ -130,8 +186,8 @@ Return type of element of the array that will represent the result of function `
     @inline eti_to_one(eti) = one(number_eltype(eti))
     return typeof(sum(map(eti_to_one, args)))
 end
-@inline function allocate_result_type(::AbstractManifold, f::TF, args::Tuple{}) where {TF}
-    return Float64
+@inline function allocate_result_type(M::AbstractManifold, f::TF, args::Tuple{}) where {TF}
+    return allocation_promotion_function(M, f, ())(Float64)
 end
 
 """
@@ -351,6 +407,21 @@ function at the level, where also information from `p` and `M` can be accessed.
 """
 copyto!(::AbstractManifold, Y, p, X) = copyto!(Y, X)
 
+"""
+    default_type(M::AbstractManifold)
+
+Get the default type of points on manifold `M`.
+"""
+default_type(M::AbstractManifold) = typeof(allocate_on(M))
+"""
+    default_type(M::AbstractManifold, ft::FiberType)
+
+Get the default type of points from the fiber `ft` of the fiber bundle based on manifold `M`.
+For example, call `default_type(MyManifold(), TangentSpaceType())` to get the default type
+of a tangent vector.
+"""
+default_type(M::AbstractManifold, ft::FiberType) = typeof(allocate_on(M, ft))
+
 @doc raw"""
     distance(M::AbstractManifold, p, q)
 
@@ -477,7 +548,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 +558,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))
@@ -1115,6 +1186,9 @@ function __init__()
         @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin
             include("../ext/ManifoldsBasePlotsExt.jl")
         end
+        @require Quaternions = "94ee1d12-ae83-5a48-8b1c-48b8ff168ae0" begin
+            include("../ext/ManifoldsBaseQuaternionsExt.jl")
+        end
         @require Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" begin
             include("../ext/ManifoldsBaseStatisticsExt.jl")
         end
@@ -1219,6 +1293,7 @@ export ×,
     ℝ,
     ℂ,
     allocate,
+    allocate_on,
     angle,
     base_manifold,
     base_point,
@@ -1236,6 +1311,7 @@ export ×,
     default_approximation_method,
     default_inverse_retraction_method,
     default_retraction_method,
+    default_type,
     default_vector_transport_method,
     distance,
     exp,

src/PowerManifold.jl

@@ -172,6 +172,52 @@ function Base.:^(
     return PowerManifold(M, size...)
 end
 
+function allocate_on(
+    M::PowerManifold{
+        𝔽,
+        TM,
+        TSize,
+        <:Union{NestedPowerRepresentation,NestedReplacingPowerRepresentation},
+    },
+) where {𝔽,TM<:AbstractManifold{𝔽},TSize}
+    return [allocate_on(M.manifold) for _ in get_iterator(M)]
+end
+function allocate_on(
+    M::PowerManifold{
+        𝔽,
+        TM,
+        TSize,
+        <:Union{NestedPowerRepresentation,NestedReplacingPowerRepresentation},
+    },
+    ::Type{<:Array{U}},
+) where {𝔽,TM<:AbstractManifold{𝔽},TSize,U}
+    return [allocate_on(M.manifold, U) for _ in get_iterator(M)]
+end
+
+function allocate_on(
+    M::PowerManifold{
+        𝔽,
+        TM,
+        TSize,
+        <:Union{NestedPowerRepresentation,NestedReplacingPowerRepresentation},
+    },
+    ft::TangentSpaceType,
+) where {𝔽,TM<:AbstractManifold{𝔽},TSize}
+    return [allocate_on(M.manifold, ft) for _ in get_iterator(M)]
+end
+function allocate_on(
+    M::PowerManifold{
+        𝔽,
+        TM,
+        TSize,
+        <:Union{NestedPowerRepresentation,NestedReplacingPowerRepresentation},
+    },
+    ft::TangentSpaceType,
+    ::Type{<:Array{U}},
+) where {𝔽,TM<:AbstractManifold{𝔽},TSize,U}
+    return [allocate_on(M.manifold, ft, U) for _ in get_iterator(M)]
+end
+
 """
     _allocate_access_nested(M::PowerManifoldNested, y, i)
 
@@ -194,7 +240,7 @@ function allocate_result(M::PowerManifoldNested, f, x...)
         ]
     end
 end
-# avoid ambituities - though usually not used
+# avoid ambiguities - though usually not used
 function allocate_result(
     M::PowerManifoldNested,
     f::typeof(get_coordinates),

src/bases.jl

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