Basis number system swap (#754)

* adapt Euclidean to basis number swap * fix symmetric and unitary * fix projective space * BVP adjustments * the new BVP API fails on Julia 1.6 * add date
JuliaManifolds · Oct 4, 2024 · 029b766 · 029b766
1 parent 150e7e0
commit 029b766
Show file tree

Hide file tree

Showing 12 changed files with 118 additions and 73 deletions.
diff --git a/NEWS.md b/NEWS.md
@@ -5,12 +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.10.3] - unreleased
+## [0.10.3] - 2024-10-04
+
+### Changed
+
+* **Mildly breaking**: the number system parameter now corresponds to the coefficients standing in front of basis vectors in a linear combination instead of components of a vector. For example, `DefaultOrthonormalBasis() == DefaultOrthonormalBasis(ℝ)` of `Euclidean(3, field=ℂ)` now has 6 vectors, and `DefaultOrthonormalBasis(ℂ)` of the same manifold has 3 basis vectors.
+
+### Fixed
 
 * Fixed `solve_exp_ode` only returning the starting position ([#744](https://github.com/JuliaManifolds/Manifolds.jl/issues/744))
 * Fixed documentation of `solve_exp_ode` function signature ([#740](https://github.com/JuliaManifolds/Manifolds.jl/issues/740))  
 
-## [0.10.2] - unreleased
+## [0.10.2] - 2024-09-24
 
 ### Added
 

diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "Manifolds"
 uuid = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e"
 authors = ["Seth Axen <[email protected]>", "Mateusz Baran <[email protected]>", "Ronny Bergmann <[email protected]>", "Antoine Levitt <[email protected]>"]
-version = "0.10.2"
+version = "0.10.3"
 
 [deps]
 Einsum = "b7d42ee7-0b51-5a75-98ca-779d3107e4c0"
@@ -55,7 +55,7 @@ HybridArrays = "0.4"
 Kronecker = "0.4, 0.5"
 LinearAlgebra = "1.6"
 ManifoldDiff = "0.3.7"
-ManifoldsBase = "0.15.15"
+ManifoldsBase = "0.15.17"
 Markdown = "1.6"
 MatrixEquations = "2.2"
 NLsolve = "4"

diff --git a/ext/ManifoldsTestExt/tests_general.jl b/ext/ManifoldsTestExt/tests_general.jl
@@ -651,10 +651,7 @@ function test_manifold(
         q = pts[2]
         X = inverse_retract(M, p, q, default_inverse_retraction_method)
         Y = vee(M, p, X)
-        Test.@test length(Y) == number_of_coordinates(
-            M,
-            ManifoldsBase.VeeOrthogonalBasis(number_system(M)),
-        )
+        Test.@test length(Y) == number_of_coordinates(M, ManifoldsBase.VeeOrthogonalBasis())
         Test.@test isapprox(M, p, X, hat(M, p, Y))
         Y2 = allocate(Y)
         vee_ret = vee!(M, Y2, p, X)

diff --git a/src/manifolds/Euclidean.jl b/src/manifolds/Euclidean.jl
@@ -226,11 +226,18 @@ function get_basis_diagonalizing(
     return CachedBasis(B, DiagonalizingBasisData(B.frame_direction, eigenvalues, vecs))
 end
 
-function get_coordinates_orthonormal(::Euclidean, p, X, ::RealNumbers)
+function get_coordinates_orthonormal(::Euclidean{<:Any,ℝ}, p, X, ::RealNumbers)
+    return vec(X)
+end
+function get_coordinates_orthonormal(::Euclidean{<:Any,ℂ}, p, X, ::ComplexNumbers)
     return vec(X)
 end
 
-function get_coordinates_orthonormal!(::Euclidean, c, p, X, ::RealNumbers)
+function get_coordinates_orthonormal!(::Euclidean{<:Any,ℝ}, c, p, X, ::RealNumbers)
+    copyto!(c, vec(X))
+    return c
+end
+function get_coordinates_orthonormal!(::Euclidean{<:Any,ℂ}, c, p, X, ::ComplexNumbers)
     copyto!(c, vec(X))
     return c
 end
@@ -248,7 +255,7 @@ function get_coordinates_induced_basis!(
     return c
 end
 
-function get_coordinates_orthonormal!(M::Euclidean{<:Any,ℂ}, c, ::Any, X, ::ComplexNumbers)
+function get_coordinates_orthonormal!(M::Euclidean{<:Any,ℂ}, c, ::Any, X, ::RealNumbers)
     S = representation_size(M)
     PS = prod(S)
     c .= [reshape(real.(X), PS)..., reshape(imag(X), PS)...]
@@ -260,27 +267,31 @@ function get_coordinates_diagonalizing!(
     c,
     ::Any,
     X,
-    ::DiagonalizingOrthonormalBasis{ℂ},
+    ::DiagonalizingOrthonormalBasis{ℝ},
 )
     S = representation_size(M)
     PS = prod(S)
     c .= [reshape(real.(X), PS)..., reshape(imag(X), PS)...]
     return c
 end
 function get_coordinates_diagonalizing!(
-    M::Euclidean,
+    M::Euclidean{<:Any,𝔽},
     c,
     p,
     X,
-    ::DiagonalizingOrthonormalBasis{ℝ},
-)
+    ::DiagonalizingOrthonormalBasis{𝔽},
+) where {𝔽}
     S = representation_size(M)
     PS = prod(S)
     copyto!(c, reshape(X, PS))
     return c
 end
 
-function get_vector_orthonormal(M::Euclidean, ::Any, c, ::RealNumbers)
+function get_vector_orthonormal(M::Euclidean{<:Any,ℝ}, ::Any, c, ::RealNumbers)
+    S = representation_size(M)
+    return reshape(c, S)
+end
+function get_vector_orthonormal(M::Euclidean{<:Any,ℂ}, ::Any, c, ::ComplexNumbers)
     S = representation_size(M)
     return reshape(c, S)
 end
@@ -298,7 +309,7 @@ function get_vector_orthonormal(::Euclidean{Tuple{Int},ℝ}, ::Any, c, ::RealNum
     return c
 end
 function get_vector_orthonormal(
-    ::Euclidean{<:TypeParameter},
+    ::Euclidean{<:TypeParameter,ℝ},
     ::SArray{S},
     c,
     ::RealNumbers,
@@ -314,14 +325,6 @@ function get_vector_orthonormal(
     # probably doesn't need rewrapping in SArray
     return c
 end
-function Manifolds.get_vector_orthonormal(
-    ::Euclidean{TypeParameter{Tuple{N}}},
-    ::SizedArray{S},
-    c,
-    ::RealNumbers,
-) where {N,S}
-    return SizedArray{S}(c)
-end
 function get_vector_orthonormal(
     ::Euclidean{TypeParameter{Tuple{N}},ℝ},
     ::SizedArray{S},
@@ -343,7 +346,12 @@ function get_vector_orthonormal!(
     copyto!(Y, c)
     return Y
 end
-function get_vector_orthonormal!(M::Euclidean, Y, ::Any, c, ::RealNumbers)
+function get_vector_orthonormal!(M::Euclidean{<:Any,ℝ}, Y, ::Any, c, ::RealNumbers)
+    S = representation_size(M)
+    copyto!(Y, reshape(c, S))
+    return Y
+end
+function get_vector_orthonormal!(M::Euclidean{<:Any,ℂ}, Y, ::Any, c, ::ComplexNumbers)
     S = representation_size(M)
     copyto!(Y, reshape(c, S))
     return Y
@@ -364,7 +372,7 @@ function get_vector_induced_basis!(M::Euclidean, Y, ::Any, c, B::InducedBasis)
     copyto!(Y, reshape(c, S))
     return Y
 end
-function get_vector_orthonormal!(M::Euclidean{<:Any,ℂ}, Y, ::Any, c, ::ComplexNumbers)
+function get_vector_orthonormal!(M::Euclidean{<:Any,ℂ}, Y, ::Any, c, ::RealNumbers)
     S = representation_size(M)
     N = div(length(c), 2)
     copyto!(Y, reshape(c[1:N] + im * c[(N + 1):end], S))
@@ -375,7 +383,7 @@ function get_vector_diagonalizing!(
     Y,
     ::Any,
     c,
-    ::DiagonalizingOrthonormalBasis{ℂ},
+    ::DiagonalizingOrthonormalBasis{ℝ},
 )
     S = representation_size(M)
     N = div(length(c), 2)

diff --git a/src/manifolds/ProjectiveSpace.jl b/src/manifolds/ProjectiveSpace.jl
@@ -207,21 +207,38 @@ formula for ``Y`` is
 """
 get_coordinates(::AbstractProjectiveSpace{ℝ}, p, X, ::DefaultOrthonormalBasis)
 
-function get_coordinates_orthonormal!(
-    M::AbstractProjectiveSpace{𝔽},
-    Y,
-    p,
-    X,
-    ::RealNumbers,
-) where {𝔽}
-    n = div(manifold_dimension(M), real_dimension(𝔽))
+function _gc_impl!(c, p, X, n::Int)
     z = p[1]
     cosθ = abs(z)
     λ = nzsign(z, cosθ)
     pend, Xend = view(p, 2:(n + 1)), view(X, 2:(n + 1))
     factor = λ' * X[1] / (1 + cosθ)
-    Y .= (Xend .- pend .* factor) .* λ'
-    return Y
+    c .= (Xend .- pend .* factor) .* λ'
+    return c
+end
+function get_coordinates_orthonormal!(M::AbstractProjectiveSpace{ℝ}, c, p, X, ::RealNumbers)
+    n = manifold_dimension(M)
+    return _gc_impl!(c, p, X, n)
+end
+function get_coordinates_orthonormal!(
+    M::AbstractProjectiveSpace{ℂ},
+    c,
+    p,
+    X,
+    ::ComplexNumbers,
+)
+    n = div(manifold_dimension(M), 2)
+    return _gc_impl!(c, p, X, n)
+end
+function get_coordinates_orthonormal!(
+    M::AbstractProjectiveSpace{ℍ},
+    c,
+    p,
+    X,
+    ::QuaternionNumbers,
+)
+    n = div(manifold_dimension(M), 4)
+    return _gc_impl!(c, p, X, n)
 end
 
 @doc raw"""
@@ -242,14 +259,7 @@ Y = \left(X - q\frac{2 \left\langle q, \begin{pmatrix}0 \\ X\end{pmatrix}\right\
 """
 get_vector(::AbstractProjectiveSpace, p, X, ::DefaultOrthonormalBasis{ℝ})
 
-function get_vector_orthonormal!(
-    M::AbstractProjectiveSpace{𝔽},
-    Y,
-    p,
-    X,
-    ::RealNumbers,
-) where {𝔽}
-    n = div(manifold_dimension(M), real_dimension(𝔽))
+function _gv_impl!(Y, p, X, n::Int)
     z = p[1]
     cosθ = abs(z)
     λ = nzsign(z, cosθ)
@@ -259,6 +269,24 @@ function get_vector_orthonormal!(
     Y[2:(n + 1)] .= (X .- pend .* (pX / (1 + cosθ))) .* λ
     return Y
 end
+function get_vector_orthonormal!(M::AbstractProjectiveSpace{ℝ}, Y, p, X, ::RealNumbers)
+    n = manifold_dimension(M)
+    return _gv_impl!(Y, p, X, n)
+end
+function get_vector_orthonormal!(M::AbstractProjectiveSpace{ℂ}, Y, p, X, ::ComplexNumbers)
+    n = div(manifold_dimension(M), 2)
+    return _gv_impl!(Y, p, X, n)
+end
+function get_vector_orthonormal!(
+    M::AbstractProjectiveSpace{ℍ},
+    Y,
+    p,
+    X,
+    ::QuaternionNumbers,
+)
+    n = div(manifold_dimension(M), 4)
+    return _gv_impl!(Y, p, X, n)
+end
 
 injectivity_radius(::AbstractProjectiveSpace) = π / 2
 injectivity_radius(::AbstractProjectiveSpace, p) = π / 2

diff --git a/src/manifolds/Symmetric.jl b/src/manifolds/Symmetric.jl
@@ -103,13 +103,7 @@ function get_coordinates_orthonormal!(M::SymmetricMatrices{<:Any,ℝ}, Y, p, X,
     end
     return Y
 end
-function get_coordinates_orthonormal!(
-    M::SymmetricMatrices{<:Any,ℂ},
-    Y,
-    p,
-    X,
-    ::ComplexNumbers,
-)
+function get_coordinates_orthonormal!(M::SymmetricMatrices{<:Any,ℂ}, Y, p, X, ::RealNumbers)
     N = get_parameter(M.size)[1]
     dim = manifold_dimension(M)
     @assert size(Y) == (dim,)
@@ -150,7 +144,7 @@ function get_vector_orthonormal!(M::SymmetricMatrices{<:Any,ℝ}, Y, p, X, ::Rea
     end
     return Y
 end
-function get_vector_orthonormal!(M::SymmetricMatrices{<:Any,ℂ}, Y, p, X, ::ComplexNumbers)
+function get_vector_orthonormal!(M::SymmetricMatrices{<:Any,ℂ}, Y, p, X, ::RealNumbers)
     N = get_parameter(M.size)[1]
     dim = manifold_dimension(M)
     @assert size(X) == (dim,)

diff --git a/src/manifolds/Unitary.jl b/src/manifolds/Unitary.jl
@@ -53,7 +53,7 @@ function get_coordinates_orthonormal(
     ::UnitaryMatrices{TypeParameter{Tuple{1}},ℍ},
     p,
     X::Quaternions.Quaternion,
-    ::QuaternionNumbers,
+    ::RealNumbers,
 )
     return @SVector [X.v1, X.v2, X.v3]
 end
@@ -62,7 +62,7 @@ function get_vector_orthonormal(
     ::UnitaryMatrices{TypeParameter{Tuple{1}},ℍ},
     p::Quaternions.Quaternion,
     c,
-    ::QuaternionNumbers,
+    ::RealNumbers,
 )
     i = firstindex(c)
     return Quaternions.quat(0, c[i], c[i + 1], c[i + 2])

diff --git a/test/manifolds/embedded_torus.jl b/test/manifolds/embedded_torus.jl
@@ -47,9 +47,21 @@ using BoundaryValueDiffEq
         final_time=3.0,
     )
     pexp_3 = p_exp(3.0)
-    @test pexp_3[1] ≈ [2.701765894057119, 2.668437820810143, -1.8341712552932237]
-    @test pexp_3[2] ≈ [-0.41778834843865575, 2.935021992911625, 0.7673987137187901]
-    @test pexp_3[3] ≈ [7.661627684089519, 4.629037950515605, 3.7839234533367194]
+    @test isapprox(
+        pexp_3[1],
+        [2.701765894057119, 2.668437820810143, -1.8341712552932237];
+        atol=1e-5,
+    )
+    @test isapprox(
+        pexp_3[2],
+        [-0.41778834843865575, 2.935021992911625, 0.7673987137187901];
+        atol=1e-5,
+    )
+    @test isapprox(
+        pexp_3[3],
+        [7.661627684089519, 4.629037950515605, 3.7839234533367194];
+        atol=1e-4,
+    )
     @test_throws DomainError p_exp(10.0)
     @test_throws DomainError p_exp(-1.0)
     @test p_exp([3.0]) == [pexp_3]

diff --git a/test/manifolds/euclidean.jl b/test/manifolds/euclidean.jl
@@ -364,7 +364,7 @@ using FiniteDifferences
             M1c,
             SizedVector{3}([1.0im, 2.0, 4.0im]),
             SizedVector{3}([-1.0, -3.0, -4.0im]),
-            DefaultOrthonormalBasis(),
+            DefaultOrthonormalBasis(ℂ),
         ) == SA[-1.0, -3.0, -4.0im]
     end
 

diff --git a/test/manifolds/projective_space.jl b/test/manifolds/projective_space.jl
@@ -137,12 +137,12 @@ include("../header.jl")
         end
         types = [Vector{ComplexF64}]
         @testset "Type $T" for T in types
-            x = [0.5 + 0.5im, 0.5 + 0.5im, 0]
-            v = [0.0, 0.0, 1.0 - im]
-            y = im * exp(M, x, v)
-            w = [0.5, -0.5, 0.5im]
-            z = (sqrt(0.5) - sqrt(0.5) * im) * exp(M, x, w)
-            pts = convert.(T, [x, y, z])
+            p1 = [0.5 + 0.5im, 0.5 + 0.5im, 0]
+            X = [0.0, 0.0, 1.0 - im]
+            p2 = im * exp(M, p1, X)
+            Y = [0.5, -0.5, 0.5im]
+            p3 = (sqrt(0.5) - sqrt(0.5) * im) * exp(M, p1, Y)
+            pts = convert.(T, [p1, p2, p3])
             test_manifold(
                 M,
                 pts,
@@ -157,7 +157,7 @@ include("../header.jl")
                     SchildsLadderTransport(),
                     PoleLadderTransport(),
                 ],
-                basis_types_to_from=(DefaultOrthonormalBasis(),),
+                basis_types_to_from=(DefaultOrthonormalBasis(ℂ),),
                 test_vee_hat=false,
                 retraction_methods=[
                     ProjectionRetraction(),
@@ -254,7 +254,7 @@ include("../header.jl")
                     SchildsLadderTransport(),
                     PoleLadderTransport(),
                 ],
-                basis_types_to_from=(DefaultOrthonormalBasis(),),
+                basis_types_to_from=(DefaultOrthonormalBasis(ℍ),),
                 test_vee_hat=false,
                 retraction_methods=[
                     ProjectionRetraction(),

diff --git a/test/manifolds/symmetric.jl b/test/manifolds/symmetric.jl
@@ -77,8 +77,8 @@ include("../header.jl")
                     SchildsLadderTransport(),
                     PoleLadderTransport(),
                 ],
-                basis_types_vecs=(DefaultOrthonormalBasis(ℂ),),
-                basis_types_to_from=(DefaultOrthonormalBasis(ℂ),),
+                basis_types_vecs=(DefaultOrthonormalBasis(ℝ),),
+                basis_types_to_from=(DefaultOrthonormalBasis(ℝ),),
                 is_tangent_atol_multiplier=1,
                 test_inplace=true,
             )