From 029b766a67f2ebf1ded8458de86a4f4f513a48e2 Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Fri, 4 Oct 2024 19:45:33 +0200 Subject: [PATCH] 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 --- NEWS.md | 10 ++++- Project.toml | 4 +- ext/ManifoldsTestExt/tests_general.jl | 5 +-- src/manifolds/Euclidean.jl | 48 +++++++++++--------- src/manifolds/ProjectiveSpace.jl | 64 +++++++++++++++++++-------- src/manifolds/Symmetric.jl | 10 +---- src/manifolds/Unitary.jl | 4 +- test/manifolds/embedded_torus.jl | 18 ++++++-- test/manifolds/euclidean.jl | 2 +- test/manifolds/projective_space.jl | 16 +++---- test/manifolds/symmetric.jl | 4 +- test/manifolds/unitary_matrices.jl | 6 +-- 12 files changed, 118 insertions(+), 73 deletions(-) diff --git a/NEWS.md b/NEWS.md index cc28716c56..95c19c0631 100644 --- 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 index ba498eba26..074f5a4b82 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Manifolds" uuid = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e" authors = ["Seth Axen ", "Mateusz Baran ", "Ronny Bergmann ", "Antoine Levitt "] -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 index cc056bff49..0607bc314f 100644 --- 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 index aec3f9c25f..9b421c6e60 100644 --- 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,7 +267,7 @@ function get_coordinates_diagonalizing!( c, ::Any, X, - ::DiagonalizingOrthonormalBasis{ℂ}, + ::DiagonalizingOrthonormalBasis{ℝ}, ) S = representation_size(M) PS = prod(S) @@ -268,19 +275,23 @@ function get_coordinates_diagonalizing!( 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 index 6d3ea67d4a..3de9320596 100644 --- 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 index 4b88fe9141..b3e4dc572d 100644 --- 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 index fb2412b27e..96aa6c6af7 100644 --- 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 index b8720d10d5..a24ab83526 100644 --- 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 index 0099f030ec..44673e6f48 100644 --- 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 index 8a843d6096..a72ff1e709 100644 --- 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 index 79a164dc97..94feb505e0 100644 --- 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, ) diff --git a/test/manifolds/unitary_matrices.jl b/test/manifolds/unitary_matrices.jl index ccd5a889f3..d0fbe27a94 100644 --- a/test/manifolds/unitary_matrices.jl +++ b/test/manifolds/unitary_matrices.jl @@ -144,12 +144,12 @@ end @test is_point(M, q) @test isapprox(q, sign(pu)) - @test get_coordinates(M, p, Quaternion(0, 1, 2, 3), DefaultOrthonormalBasis(ℍ)) == + @test get_coordinates(M, p, Quaternion(0, 1, 2, 3), DefaultOrthonormalBasis(ℝ)) == SA[1, 2, 3] - @test get_vector(M, p, SA[1, 2, 3], DefaultOrthonormalBasis(ℍ)) == + @test get_vector(M, p, SA[1, 2, 3], DefaultOrthonormalBasis(ℝ)) == Quaternion(0, 1, 2, 3) - @test get_basis(M, p, DefaultOrthonormalBasis(ℍ)).data == [ + @test get_basis(M, p, DefaultOrthonormalBasis(ℝ)).data == [ Quaternion(0.0, 1.0, 0.0, 0.0), Quaternion(0.0, 0.0, 1.0, 0.0), Quaternion(0.0, 0.0, 0.0, 1.0),