From bd0c6ba0a3d9908977f41e55e0992a0061692ed6 Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Sat, 16 Nov 2019 19:01:03 +0100 Subject: [PATCH 01/12] rename Euclidean to DefaultManifold and shorten the file to not include documentation and extends the tests to also include geodesics. --- src/DefaultManifold.jl | 25 ++++++++++ src/ManifoldsBase.jl | 6 ++- src/euclidean.jl | 54 ---------------------- test/{euclidean.jl => default_manifold.jl} | 23 +++++++-- test/runtests.jl | 2 +- 5 files changed, 50 insertions(+), 60 deletions(-) create mode 100644 src/DefaultManifold.jl delete mode 100644 src/euclidean.jl rename test/{euclidean.jl => default_manifold.jl} (82%) diff --git a/src/DefaultManifold.jl b/src/DefaultManifold.jl new file mode 100644 index 00000000..066db98e --- /dev/null +++ b/src/DefaultManifold.jl @@ -0,0 +1,25 @@ +struct DefaultManifold{T<:Tuple} <: Manifold where {T} end +DefaultManifold(n::Vararg{Int,N}) where N = DefaultManifold{Tuple{n...}}() +@generated representation_size(::DefaultManifold{T}) where {T} = Tuple(T.parameters...) +@generated manifold_dimension(::DefaultManifold{T}) where {T} = *(T.parameters...) +@inline inner(::DefaultManifold, x, v, w) = dot(v, w) +distance(::DefaultManifold, x, y) = norm(x-y) +norm(::DefaultManifold, x, v) = norm(v) +exp!(M::DefaultManifold, y, x, v) = (y .= x .+ v) +log!(M::DefaultManifold, v, x, y) = (v .= y .- x) +function zero_tangent_vector!(M::DefaultManifold, v, x) + fill!(v, 0) + return v +end +function project_point!(M::DefaultManifold, y, x) + y .= x + return y +end +function project_tangent!(M::DefaultManifold, w, x, v) + w .= v + return w +end +function vector_transport_to!(M::DefaultManifold, vto, x, v, y, ::ParallelTransport) + vto .= v + return vto +end \ No newline at end of file diff --git a/src/ManifoldsBase.jl b/src/ManifoldsBase.jl index 2a093b76..2c9aa6c7 100644 --- a/src/ManifoldsBase.jl +++ b/src/ManifoldsBase.jl @@ -604,10 +604,12 @@ function check_tangent_vector(M::Manifold, x, v; kwargs...) end -include("euclidean.jl") +include("DefaultManifold.jl") export Manifold, - Euclidean + MPoint, + TVector, + CoTVector export base_manifold, check_manifold_point, diff --git a/src/euclidean.jl b/src/euclidean.jl deleted file mode 100644 index c69cdc76..00000000 --- a/src/euclidean.jl +++ /dev/null @@ -1,54 +0,0 @@ -@doc doc""" - Euclidean{T<:Tuple} <: Manifold - -Euclidean vector space $\mathbb R^n$. - -# Constructor - - Euclidean(n) - -generates the $n$-dimensional vector space $\mathbb R^n$. - - Euclidean(m, n) - -generates the $mn$-dimensional vector space $\mathbb R^{m \times n}$, whose -elements are interpreted as $m \times n$ matrices. -""" -struct Euclidean{T<:Tuple} <: Manifold where {T} end - -Euclidean(n::Int) = Euclidean{Tuple{n}}() -Euclidean(m::Int, n::Int) = Euclidean{Tuple{m,n}}() - -function representation_size(::Euclidean{Tuple{n}}) where {n} - return (n,) -end - -function representation_size(::Euclidean{Tuple{m,n}}) where {m,n} - return (m,n) -end - -@generated manifold_dimension(::Euclidean{T}) where {T} = *(T.parameters...) - -@inline inner(::Euclidean, x, v, w) = dot(v, w) - -distance(::Euclidean, x, y) = norm(x-y) -norm(::Euclidean, x, v) = norm(v) - -exp!(M::Euclidean, y, x, v) = (y .= x .+ v) - -log!(M::Euclidean, v, x, y) = (v .= y .- x) - -function zero_tangent_vector!(M::Euclidean, v, x) - fill!(v, 0) - return v -end - -function project_point!(M::Euclidean, y, x) - copyto!(y, x) - return y -end - -function project_tangent!(M::Euclidean, w, x, v) - w .= v - return w -end diff --git a/test/euclidean.jl b/test/default_manifold.jl similarity index 82% rename from test/euclidean.jl rename to test/default_manifold.jl index 2ef0f486..9fbc54ef 100644 --- a/test/euclidean.jl +++ b/test/default_manifold.jl @@ -7,8 +7,8 @@ using ReverseDiff using StaticArrays using Test -@testset "Euclidean" begin - M = ManifoldsBase.Euclidean(3) +@testset "Testting Default (Euclidean)" begin + M = ManifoldsBase.DefaultManifold(3) types = [Vector{Float64}, SizedVector{3, Float64}, MVector{3, Float64}, @@ -34,7 +34,6 @@ using Test pts = [convert(T, [1.0, 0.0, 0.0]), convert(T, [0.0, 1.0, 0.0]), convert(T, [0.0, 0.0, 1.0])] - @test injectivity_radius(M, pts[1]) == Inf @test injectivity_radius(M, pts[1], rm) == Inf @@ -74,6 +73,24 @@ using Test @test distance(M, pts[1], pts[2]) ≈ norm(M, pts[1], tv1) + @testset "Geodesic Interface Test" begin + @test isapprox(M, geodesic(M, pts[1], tv1)(0.), pts[1]) + @test isapprox(M, geodesic(M, pts[1], tv1)(1.), pts[2]) + @test isapprox(M, geodesic(M, pts[1],tv1, 1.), pts[2]) + @test isapprox(M, geodesic(M, pts[1],tv1, 1. /2), (pts[1]+pts[2])/2) + @test isapprox(M, shortest_geodesic(M, pts[1], pts[2])(0.), pts[1]) + @test isapprox(M, shortest_geodesic(M, pts[1], pts[2])(1.), pts[2]) + @test isapprox(M, shortest_geodesic(M, pts[1], pts[2], 0.), pts[1]) + @test isapprox(M, shortest_geodesic(M, pts[1], pts[2], 1.), pts[2]) + @test all( + isapprox.(Ref(M), geodesic(M, pts[1], tv1, [0., 1. /2, 1.]), + [pts[1], (pts[1]+pts[2])/2, pts[2]] ) + ) + @test all( + isapprox.(Ref(M), shortest_geodesic(M, pts[1], pts[2], [0., 1. /2, 1.]), + [pts[1], (pts[1]+pts[2])/2, pts[2]] ) + ) + end @testset "basic linear algebra in tangent space" begin @test isapprox(M, pts[1], 0*tv1, zero_tangent_vector(M, pts[1]); atol = eps(eltype(pts[1]))) diff --git a/test/runtests.jl b/test/runtests.jl index 0de11238..33a646cc 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,3 +1,3 @@ include("empty_manifold.jl") -include("euclidean.jl") +include("default_manifold.jl") From ff99d9403a84706ae9dc146991d12b25f9cac2d8 Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Sat, 16 Nov 2019 19:01:35 +0100 Subject: [PATCH 02/12] Extends the empty_tests and adds more details on vector_transports. --- README.md | 6 ++- src/ManifoldsBase.jl | 89 ++++++++++++++++++++++++++++++------------ test/empty_manifold.jl | 63 ++++++++++++++++-------------- 3 files changed, 104 insertions(+), 54 deletions(-) diff --git a/README.md b/README.md index 57b9157d..5a8abc56 100644 --- a/README.md +++ b/README.md @@ -4,4 +4,8 @@ Basic interface for manifolds in Julia. -The main part of the project is [`Manifolds.jl`](https://github.com/JuliaNLSolvers/Manifolds.jl). +This interface includes an easy `DefaultManifold`, which is a reduced version +of the `Eucliean` manifold from [`Manifolds.jl`](https://github.com/JuliaNLSolvers/Manifolds.jl) such that the interface functions can be tested. + +The project [`Manifolds.jl`](https://github.com/JuliaNLSolvers/Manifolds.jl) +is based on this interface and provides a variety of manifolds. diff --git a/src/ManifoldsBase.jl b/src/ManifoldsBase.jl index 2c9aa6c7..958eaf52 100644 --- a/src/ManifoldsBase.jl +++ b/src/ManifoldsBase.jl @@ -384,69 +384,110 @@ function shortest_geodesic(M::Manifold, x, y, T::AbstractVector) return geodesic(M, x, log(M, x, y), T) end +abstract type AbstractVectorTransportMethod end + +""" + ParallelTransport <: AbstractVectorTransportMethod + +Specify to use parallel transport as vector transport method within +[`vector_transport_to`](@ref), [`vector_transport_direction`](@ref) or +[`vector_transport_along`](@ref). +""" +struct ParallelTransport <: AbstractVectorTransportMethod end + +""" + ProjectTangent <: AbstractVectorTransportMethod + +Specify to use projection onto tangent space as vector transport method within +[`vector_transport_to`](@ref), [`vector_transport_direction`](@ref) or +[`vector_transport_along`](@ref). See [`project_tangent`](@ref) for details. +""" +struct ProjectTangent <: AbstractVectorTransportMethod end """ - vector_transport_to!(M::Manifold, vto, x, v, y) + vector_transport_to!(M::Manifold, vto, x, v, y, m::AbstractVectorTransportMethod=ParallelTransport()) Vector transport of vector `v` at point `x` to point `y`. The result is saved -to `vto`. By default, [`project_tangent!`](@ref) is used but this may change in -the future. +to `vto`. By default, the method `m` is [`ParallelTransport`](@ref). """ -vector_transport_to!(M::Manifold, vto, x, v, y) = project_tangent!(M, vto, y, v) +vector_transport_to!(M::Manifold, vto, x, v, y) = vector_transport_to!(M,vto,x,v,y,ParallelTransport()) + +""" + vector_transport_to!(M::Manifold, vto, x, v, y, ProjectTangent()) + +Implements a default projection based vector transport, that projects a tangent +vector `v` at `x` on a [`Manifold`](@ref) `M` onto the tangent space at `y` by +interperting `v` as an element of the embedding and projecting back. +""" +vector_transport_to!(M::Manifold, vto, x, v, y, m::ProjectTangent) = project_tangent!(M, vto, y, v) + +function vector_transport_to!(M::Manifold, vto, x, v, y, m::AbstractVectorTransportMethod) + error("vector transport from a point of type $(typeof(x)) to a type $(typeof(y)) on a $(typeof(M)) for a vector of type $(v) and the $(typeof(m)) not yet implemented.") +end + """ - vector_transport_to(M::Manifold, x, v, y) + vector_transport_to(M::Manifold, x, v, y, m::AbstractVectorTransportMethod=ParallelTransport()) -Vector transport of vector `v` at point `x` to point `y`. +Vector transport of vector `v` at point `x` to point `y` using the method `m`, +which defaults to [`ParallelTransport`](@ref). """ -function vector_transport_to(M::Manifold, x, v, y) +vector_transport_to(M::Manifold, x, v, y) = vector_transport_to(M,x,v,y,ParallelTransport()) +function vector_transport_to(M::Manifold, x, v, y, m::AbstractVectorTransportMethod) vto = similar_result(M, vector_transport_to, v, x, y) - vector_transport_to!(M, vto, x, v, y) + vector_transport_to!(M, vto, x, v, y, m) return vto end """ - vector_transport_direction!(M::Manifold, vto, x, v, vdir) + vector_transport_direction!(M::Manifold, vto, x, v, vdir, m=::ParallelTransport]) Vector transport of vector `v` at point `x` in the direction indicated by the tangent vector `vdir` at point `x`. The result is saved to `vto`. -By default, `exp` and `vector_transport_to!` are used. +By default, `exp` and `vector_transport_to!` are used with the method `m` which +defaults to [`ParallelTransport`](@ref). """ -function vector_transport_direction!(M::Manifold, vto, x, v, vdir) +vector_transport_direction!(M::Manifold, vto, x, v, vdir) = vector_transport_direction!(M,vto,x,v,vdir,ParallelTransport()) +function vector_transport_direction!(M::Manifold, vto, x, v, vdir,m::AbstractVectorTransportMethod) y = exp(M, x, vdir) - return vector_transport_to!(M, vto, x, v, y) + return vector_transport_to!(M, vto, x, v, y, m) end """ - vector_transport_direction(M::Manifold, x, v, vdir) + vector_transport_direction(M::Manifold, x, v, vdir[, m=::ParallelTransport]) Vector transport of vector `v` at point `x` in the direction indicated -by the tangent vector `vdir` at point `x`. +by the tangent vector `vdir` at point `x` using the method `m`, which defaults to [`ParallelTransport`](@ref). """ -function vector_transport_direction(M::Manifold, x, v, vdir) +vector_transport_direction(M::Manifold, x, v, vdir) = vector_transport_direction(M,x,v,vdir,ParallelTransport()) +function vector_transport_direction(M::Manifold, x, v, vdir, m::AbstractVectorTransportMethod) vto = similar_result(M, vector_transport_direction, v, x, vdir) - vector_transport_direction!(M, vto, x, v, vdir) + vector_transport_direction!(M, vto, x, v, vdir,m) return vto end """ - vector_transport_along!(M::Manifold, vto, x, v, c) + vector_transport_along!(M::Manifold, vto, x, v, c[,m=::ParallelTransport]) Vector transport of vector `v` at point `x` along the curve `c` such that -`c(0)` is equal to `x` to point `c(1)`. The result is saved to `vto`. +`c(0)` is equal to `x` to point `c(1)` using the method `m`, which defaults to +[`ParallelTransport`](@ref). The result is saved to `vto`. """ -function vector_transport_along!(M::Manifold, vto, x, v, c) - error("vector_transport_along! not implemented for manifold $(typeof(M)), vector $(typeof(vto)), point $(typeof(x)), vector $(typeof(v)) and curve $(typeof(c)).") +vector_transport_along!(M::Manifold, vto, x, v, c) = vector_transport_along!(M, vto, x, v, c, ParallelTransport()) +function vector_transport_along!(M::Manifold, vto, x, v, c, m::AbstractVectorTransportMethod) + error("vector_transport_along! not implemented for manifold $(typeof(M)), vector $(typeof(vto)), point $(typeof(x)), vector $(typeof(v)) along curve $(typeof(c)) with method $(typeof(m)).") end """ - vector_transport_along(M::Manifold, x, v, c) + vector_transport_along(M::Manifold, x, v, c[,m]) Vector transport of vector `v` at point `x` along the curve `c` such that `c(0)` is equal to `x` to point `c(1)`. +The default method `m` used is [`ParallelTransport`](@ref). """ -function vector_transport_along(M::Manifold, x, v, c) - vto = similar_result(M, vector_transport_along, x, v) - vector_transport_along!(M, vto, x, v, c) +vector_transport_along(M::Manifold, x, v, c) = vector_transport_along(M,x,v,c,ParallelTransport()) +function vector_transport_along(M::Manifold, x, v, c, m::AbstractVectorTransportMethod) + vto = similar_result(M, vector_transport_along, v, x) + vector_transport_along!(M, vto, x, v, c, m) return vto end diff --git a/test/empty_manifold.jl b/test/empty_manifold.jl index 827d192a..3446ba3a 100644 --- a/test/empty_manifold.jl +++ b/test/empty_manifold.jl @@ -1,11 +1,16 @@ using ManifoldsBase using Test - +import Base: * struct NonManifold <: Manifold end - +struct NonMPoint <: MPoint end +struct NonTVector <: TVector end +struct NonCoTVector <: CoTVector end +*(t::Float64,v::NonTVector) = v @testset "Manifold with empty implementation" begin m = NonManifold() + p = NonMPoint() + v = NonTVector() @test_throws ErrorException ManifoldsBase.representation_size(m) @@ -23,42 +28,42 @@ struct NonManifold <: Manifold end exp_retr = ManifoldsBase.ExponentialRetraction() - @test_throws ErrorException retract!(m, [0], [0], [0]) - @test_throws ErrorException retract!(m, [0], [0], [0], exp_retr) - @test_throws ErrorException retract!(m, [0], [0], [0], 0.0) - @test_throws ErrorException retract!(m, [0], [0], [0], 0.0, exp_retr) - @test_throws ErrorException retract(m, [0], [0]) - @test_throws ErrorException retract(m, [0], [0], exp_retr) - @test_throws ErrorException retract(m, [0], [0], 0.0) - @test_throws ErrorException retract(m, [0], [0], 0.0, exp_retr) + @test_throws ErrorException retract!(m, p, p, v) + @test_throws ErrorException retract!(m, p, p, v, exp_retr) + @test_throws ErrorException retract!(m, p, p, [0.0], 0.) + @test_throws ErrorException retract!(m, p, p, [0.0], 0., exp_retr) + @test_throws ErrorException retract(m, [0.0], [0.0]) + @test_throws ErrorException retract(m, [0.0], [0.0], exp_retr) + @test_throws ErrorException retract(m, [0.0], [0.0], 0.0) + @test_throws ErrorException retract(m, [0.0], [0.0], 0.0, exp_retr) log_invretr = ManifoldsBase.LogarithmicInverseRetraction() - @test_throws ErrorException inverse_retract!(m, [0], [0], [0]) - @test_throws ErrorException inverse_retract!(m, [0], [0], [0], log_invretr) - @test_throws ErrorException inverse_retract(m, [0], [0]) - @test_throws ErrorException inverse_retract(m, [0], [0], log_invretr) + @test_throws ErrorException inverse_retract!(m, p, p, p) + @test_throws ErrorException inverse_retract!(m, p, p, p, log_invretr) + @test_throws ErrorException inverse_retract(m, [0.0], [0.0]) + @test_throws ErrorException inverse_retract(m, [0.0], [0.0], log_invretr) - @test_throws ErrorException project_point!(m, [0], [0]) + @test_throws ErrorException project_point!(m, p, [0]) @test_throws ErrorException project_point(m, [0]) - @test_throws ErrorException project_tangent!(m, [0], [0], [0]) - @test_throws ErrorException project_tangent(m, [0], [0]) + @test_throws ErrorException project_tangent!(m, v, p, [0.0]) + @test_throws ErrorException project_tangent(m, [0.0], [0.0]) - @test_throws ErrorException inner(m, [0], [0], [0]) - @test_throws ErrorException norm(m, [0], [0]) - @test_throws ErrorException angle(m, [0], [0], [0]) + @test_throws ErrorException inner(m, p, v, v) + @test_throws ErrorException norm(m, p, v) + @test_throws ErrorException angle(m, p, v, v) - @test_throws ErrorException distance(m, [0], [0]) + @test_throws ErrorException distance(m, [0.0], [0.0]) - @test_throws ErrorException exp!(m, [0], [0], [0]) - @test_throws ErrorException exp!(m, [0], [0], [0], 0.0) - @test_throws ErrorException exp(m, [0], [0]) - @test_throws ErrorException exp(m, [0], [0], 0.0) - @test_throws ErrorException exp(m, [0], [0], [0]) + @test_throws ErrorException exp!(m, p, p, v) + @test_throws ErrorException exp!(m, p, p, v, 0.0) + @test_throws ErrorException exp(m, [0.0], [0.0]) + @test_throws ErrorException exp(m, [0.0], [0.0], 0.0) + @test_throws ErrorException exp(m, [0.0], [0.0], [0]) - @test_throws ErrorException log!(m, [0], [0], [0]) - @test_throws ErrorException log(m, [0], [0]) + @test_throws ErrorException log!(m, v, p, p) + @test_throws ErrorException log(m, [0.0], [0.0]) @test_throws ErrorException vector_transport_to!(m, [0], [0], [0], [0]) @test_throws ErrorException vector_transport_to(m, [0], [0], [0]) @@ -75,7 +80,7 @@ struct NonManifold <: Manifold end @test_throws ErrorException zero_tangent_vector!(m, [0], [0]) @test_throws ErrorException zero_tangent_vector(m, [0]) - + @test manifold_point_error(m, [0]) === nothing @test is_manifold_point(m, [0]) @test check_manifold_point(m, [0]) == nothing From cc3ce19fd6147fff4d95be42d6c043d971ffdf79 Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Sat, 16 Nov 2019 21:59:44 +0100 Subject: [PATCH 03/12] adds the old tests again. --- test/empty_manifold.jl | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/test/empty_manifold.jl b/test/empty_manifold.jl index 3446ba3a..b5c960f4 100644 --- a/test/empty_manifold.jl +++ b/test/empty_manifold.jl @@ -32,6 +32,14 @@ struct NonCoTVector <: CoTVector end @test_throws ErrorException retract!(m, p, p, v, exp_retr) @test_throws ErrorException retract!(m, p, p, [0.0], 0.) @test_throws ErrorException retract!(m, p, p, [0.0], 0., exp_retr) + @test_throws ErrorException retract!(m, [0], [0], [0]) + @test_throws ErrorException retract!(m, [0], [0], [0], exp_retr) + @test_throws ErrorException retract!(m, [0], [0], [0], 0.0) + @test_throws ErrorException retract!(m, [0], [0], [0], 0.0, exp_retr) + @test_throws ErrorException retract(m, [0], [0]) + @test_throws ErrorException retract(m, [0], [0], exp_retr) + @test_throws ErrorException retract(m, [0], [0], 0.0) + @test_throws ErrorException retract(m, [0], [0], 0.0, exp_retr) @test_throws ErrorException retract(m, [0.0], [0.0]) @test_throws ErrorException retract(m, [0.0], [0.0], exp_retr) @test_throws ErrorException retract(m, [0.0], [0.0], 0.0) @@ -41,28 +49,45 @@ struct NonCoTVector <: CoTVector end @test_throws ErrorException inverse_retract!(m, p, p, p) @test_throws ErrorException inverse_retract!(m, p, p, p, log_invretr) + @test_throws ErrorException inverse_retract!(m, [0], [0], [0]) + @test_throws ErrorException inverse_retract!(m, [0], [0], [0], log_invretr) + @test_throws ErrorException inverse_retract(m, [0], [0]) + @test_throws ErrorException inverse_retract(m, [0], [0], log_invretr) + @test_throws ErrorException inverse_retract(m, [0.0], [0.0]) @test_throws ErrorException inverse_retract(m, [0.0], [0.0], log_invretr) @test_throws ErrorException project_point!(m, p, [0]) + @test_throws ErrorException project_point!(m, [0], [0]) @test_throws ErrorException project_point(m, [0]) @test_throws ErrorException project_tangent!(m, v, p, [0.0]) + @test_throws ErrorException project_tangent!(m, [0], [0], [0]) + @test_throws ErrorException project_tangent(m, [0], [0]) @test_throws ErrorException project_tangent(m, [0.0], [0.0]) @test_throws ErrorException inner(m, p, v, v) + @test_throws ErrorException inner(m, [0], [0], [0]) @test_throws ErrorException norm(m, p, v) + @test_throws ErrorException norm(m, [0], [0]) @test_throws ErrorException angle(m, p, v, v) + @test_throws ErrorException angle(m, [0], [0], [0]) @test_throws ErrorException distance(m, [0.0], [0.0]) @test_throws ErrorException exp!(m, p, p, v) @test_throws ErrorException exp!(m, p, p, v, 0.0) + @test_throws ErrorException exp!(m, [0], [0], [0]) + @test_throws ErrorException exp!(m, [0], [0], [0], 0.0) + @test_throws ErrorException exp(m, [0], [0]) + @test_throws ErrorException exp(m, [0], [0], 0.0) + @test_throws ErrorException exp(m, [0], [0], [0]) @test_throws ErrorException exp(m, [0.0], [0.0]) @test_throws ErrorException exp(m, [0.0], [0.0], 0.0) @test_throws ErrorException exp(m, [0.0], [0.0], [0]) @test_throws ErrorException log!(m, v, p, p) + @test_throws ErrorException log!(m, [0], [0], [0]) @test_throws ErrorException log(m, [0.0], [0.0]) @test_throws ErrorException vector_transport_to!(m, [0], [0], [0], [0]) From 277f01792bf20429130feca230513914b005c9b5 Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Sat, 16 Nov 2019 22:02:41 +0100 Subject: [PATCH 04/12] Adds a Docstring to the `DefaultManifold`, adds Types for semantic interpretation of points, tangent and cotangent vectors. --- src/DefaultManifold.jl | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/src/DefaultManifold.jl b/src/DefaultManifold.jl index 066db98e..dc56abc1 100644 --- a/src/DefaultManifold.jl +++ b/src/DefaultManifold.jl @@ -1,4 +1,33 @@ +""" + DefaultManifold <: Manifold + +This default manifold illustrates the main features of the interface and provides +a good starting point for an own manifold. It is a simplified/shortened variant +of `Euclidean` from `Manifolds.jl`. + +this manifold further illustrates, how to type your manifold points +and tangent vectors. Note that the interface does not require this, but it +might be handy in debugging and educative situations to verify correctness of +involved variabes. +""" struct DefaultManifold{T<:Tuple} <: Manifold where {T} end + +struct DefaultMPoint{T} <: MPoint where {T} + data::Array{T} +end +struct DefaultTVector{T} <: TVector where {T} + data::Array{T} +end +struct DefaultCoTVector{T} <: CoTVector where {T} + data::Array{T} +end +for T in [DefaultMPoint,DefaultTVector,DefaultCoTVector] + Base.eltype(p::T) = eltype(p.data) + Base.similar(p::T, ::Type{EltypeNew}) where EltypeNew = T(similar(p.data, EltypeNew)) + Base.similar(p::T) = T(similar(p.data)) + Base.copy(p::T) = T(copy(p.data)) +end + DefaultManifold(n::Vararg{Int,N}) where N = DefaultManifold{Tuple{n...}}() @generated representation_size(::DefaultManifold{T}) where {T} = Tuple(T.parameters...) @generated manifold_dimension(::DefaultManifold{T}) where {T} = *(T.parameters...) From d0e6b2cd20ed9f51c3119e19d8a1980ed7812b09 Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Sat, 16 Nov 2019 22:03:05 +0100 Subject: [PATCH 05/12] rename ProjectTangent to be consistent with ParallelTransport. --- src/ManifoldsBase.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/ManifoldsBase.jl b/src/ManifoldsBase.jl index 958eaf52..ffd92ed9 100644 --- a/src/ManifoldsBase.jl +++ b/src/ManifoldsBase.jl @@ -396,13 +396,13 @@ Specify to use parallel transport as vector transport method within struct ParallelTransport <: AbstractVectorTransportMethod end """ - ProjectTangent <: AbstractVectorTransportMethod + ProjectionTransport <: AbstractVectorTransportMethod Specify to use projection onto tangent space as vector transport method within [`vector_transport_to`](@ref), [`vector_transport_direction`](@ref) or [`vector_transport_along`](@ref). See [`project_tangent`](@ref) for details. """ -struct ProjectTangent <: AbstractVectorTransportMethod end +struct ProjectionTransport <: AbstractVectorTransportMethod end """ vector_transport_to!(M::Manifold, vto, x, v, y, m::AbstractVectorTransportMethod=ParallelTransport()) @@ -412,13 +412,13 @@ to `vto`. By default, the method `m` is [`ParallelTransport`](@ref). vector_transport_to!(M::Manifold, vto, x, v, y) = vector_transport_to!(M,vto,x,v,y,ParallelTransport()) """ - vector_transport_to!(M::Manifold, vto, x, v, y, ProjectTangent()) + vector_transport_to!(M::Manifold, vto, x, v, y, ProjectionTransport()) Implements a default projection based vector transport, that projects a tangent vector `v` at `x` on a [`Manifold`](@ref) `M` onto the tangent space at `y` by interperting `v` as an element of the embedding and projecting back. """ -vector_transport_to!(M::Manifold, vto, x, v, y, m::ProjectTangent) = project_tangent!(M, vto, y, v) +vector_transport_to!(M::Manifold, vto, x, v, y, m::ProjectionTransport) = project_tangent!(M, vto, y, v) function vector_transport_to!(M::Manifold, vto, x, v, y, m::AbstractVectorTransportMethod) error("vector transport from a point of type $(typeof(x)) to a type $(typeof(y)) on a $(typeof(M)) for a vector of type $(v) and the $(typeof(m)) not yet implemented.") From 8b8915d640d5fdce4bafb5798d3be342f11d915c Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Sat, 16 Nov 2019 22:14:42 +0100 Subject: [PATCH 06/12] increase code coverage. --- src/ManifoldsBase.jl | 4 +++- test/empty_manifold.jl | 3 +++ 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/src/ManifoldsBase.jl b/src/ManifoldsBase.jl index ffd92ed9..3c5360f3 100644 --- a/src/ManifoldsBase.jl +++ b/src/ManifoldsBase.jl @@ -650,7 +650,9 @@ include("DefaultManifold.jl") export Manifold, MPoint, TVector, - CoTVector + CoTVector, + ParallelTransport, + ProjectionTransport export base_manifold, check_manifold_point, diff --git a/test/empty_manifold.jl b/test/empty_manifold.jl index b5c960f4..e66dc7c7 100644 --- a/test/empty_manifold.jl +++ b/test/empty_manifold.jl @@ -92,6 +92,7 @@ struct NonCoTVector <: CoTVector end @test_throws ErrorException vector_transport_to!(m, [0], [0], [0], [0]) @test_throws ErrorException vector_transport_to(m, [0], [0], [0]) + @test_throws ErrorException vector_transport_to!(m, [0], [0], [0], ProjectionTransport()) @test_throws ErrorException vector_transport_direction!(m, [0], [0], [0], [0]) @test_throws ErrorException vector_transport_direction(m, [0], [0], [0]) @@ -107,10 +108,12 @@ struct NonCoTVector <: CoTVector end @test_throws ErrorException zero_tangent_vector(m, [0]) @test manifold_point_error(m, [0]) === nothing + @test_throws ErrorException manifold_point_error(m,p) @test is_manifold_point(m, [0]) @test check_manifold_point(m, [0]) == nothing @test tangent_vector_error(m, [0], [0]) === nothing + @test_throws ErrorException tangent_vector_error(m,p,v) @test is_tangent_vector(m, [0], [0]) @test check_tangent_vector(m, [0], [0]) == nothing From 510431d4981d54d86dfd849cbe3f3e83351ca99f Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Sun, 17 Nov 2019 07:46:35 +0100 Subject: [PATCH 07/12] add more general Default data types, update Readme. --- README.md | 6 ++++-- src/DefaultManifold.jl | 28 ++++++++-------------------- test/default_manifold.jl | 19 ++++++++++++------- 3 files changed, 24 insertions(+), 29 deletions(-) diff --git a/README.md b/README.md index 5a8abc56..11ebfa57 100644 --- a/README.md +++ b/README.md @@ -4,8 +4,10 @@ Basic interface for manifolds in Julia. -This interface includes an easy `DefaultManifold`, which is a reduced version -of the `Eucliean` manifold from [`Manifolds.jl`](https://github.com/JuliaNLSolvers/Manifolds.jl) such that the interface functions can be tested. +This interface includes a simple `DefaultManifold`, which is a reduced version +of the [`Eucliean`](https://github.com/JuliaNLSolvers/Manifolds.jl/blob/master/src/Euclidean.jl) +manifold from [`Manifolds.jl`](https://github.com/JuliaNLSolvers/Manifolds.jl) +such that the interface functions can be tested. The project [`Manifolds.jl`](https://github.com/JuliaNLSolvers/Manifolds.jl) is based on this interface and provides a variety of manifolds. diff --git a/src/DefaultManifold.jl b/src/DefaultManifold.jl index dc56abc1..c5c681b2 100644 --- a/src/DefaultManifold.jl +++ b/src/DefaultManifold.jl @@ -11,15 +11,16 @@ might be handy in debugging and educative situations to verify correctness of involved variabes. """ struct DefaultManifold{T<:Tuple} <: Manifold where {T} end +DefaultManifold(n::Vararg{Int,N}) where N = DefaultManifold{Tuple{n...}}() struct DefaultMPoint{T} <: MPoint where {T} - data::Array{T} + data::T end struct DefaultTVector{T} <: TVector where {T} - data::Array{T} + data::T end struct DefaultCoTVector{T} <: CoTVector where {T} - data::Array{T} + data::T end for T in [DefaultMPoint,DefaultTVector,DefaultCoTVector] Base.eltype(p::T) = eltype(p.data) @@ -28,7 +29,6 @@ for T in [DefaultMPoint,DefaultTVector,DefaultCoTVector] Base.copy(p::T) = T(copy(p.data)) end -DefaultManifold(n::Vararg{Int,N}) where N = DefaultManifold{Tuple{n...}}() @generated representation_size(::DefaultManifold{T}) where {T} = Tuple(T.parameters...) @generated manifold_dimension(::DefaultManifold{T}) where {T} = *(T.parameters...) @inline inner(::DefaultManifold, x, v, w) = dot(v, w) @@ -36,19 +36,7 @@ distance(::DefaultManifold, x, y) = norm(x-y) norm(::DefaultManifold, x, v) = norm(v) exp!(M::DefaultManifold, y, x, v) = (y .= x .+ v) log!(M::DefaultManifold, v, x, y) = (v .= y .- x) -function zero_tangent_vector!(M::DefaultManifold, v, x) - fill!(v, 0) - return v -end -function project_point!(M::DefaultManifold, y, x) - y .= x - return y -end -function project_tangent!(M::DefaultManifold, w, x, v) - w .= v - return w -end -function vector_transport_to!(M::DefaultManifold, vto, x, v, y, ::ParallelTransport) - vto .= v - return vto -end \ No newline at end of file +zero_tangent_vector!(M::DefaultManifold, v, x) = fill!(v, 0) +project_point!(M::DefaultManifold, y, x) = (y .= x) +project_tangent!(M::DefaultManifold, w, x, v) = (w .= v) +vector_transport_to!(M::DefaultManifold, vto, x, v, y, ::ParallelTransport) = (vto .= v) \ No newline at end of file diff --git a/test/default_manifold.jl b/test/default_manifold.jl index 9fbc54ef..6a4a56ae 100644 --- a/test/default_manifold.jl +++ b/test/default_manifold.jl @@ -31,9 +31,8 @@ using Test for T in types @testset "Type $T" begin - pts = [convert(T, [1.0, 0.0, 0.0]), - convert(T, [0.0, 1.0, 0.0]), - convert(T, [0.0, 0.0, 1.0])] + pts = convert.(Ref(T), [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]) + @test injectivity_radius(M, pts[1]) == Inf @test injectivity_radius(M, pts[1], rm) == Inf @@ -73,7 +72,7 @@ using Test @test distance(M, pts[1], pts[2]) ≈ norm(M, pts[1], tv1) - @testset "Geodesic Interface Test" begin + @testset "Geodesic interface test" begin @test isapprox(M, geodesic(M, pts[1], tv1)(0.), pts[1]) @test isapprox(M, geodesic(M, pts[1], tv1)(1.), pts[2]) @test isapprox(M, geodesic(M, pts[1],tv1, 1.), pts[2]) @@ -126,10 +125,12 @@ using Test v1 = log(M, pts[1], pts[2]) v2 = log(M, pts[1], pts[3]) v1t1 = vector_transport_to(M, pts[1], v1, pts[3]) - v1t2 = vector_transport_direction(M, pts[1], v1, v2) + v1t2 = zero(v1t1) + vector_transport_to!(M, v1t2, pts[1], v1, v2, ProjectionTransport()) + v1t3 = vector_transport_direction(M, pts[1], v1, v2) @test is_tangent_vector(M, pts[3], v1t1) - @test is_tangent_vector(M, pts[3], v1t2) - @test isapprox(M, pts[3], v1t1, v1t2) + @test is_tangent_vector(M, pts[3], v1t3) + @test isapprox(M, pts[3], v1t1, v1t3) end @testset "ForwardDiff support" begin @@ -157,6 +158,10 @@ using Test @test ReverseDiff.gradient(retract_f, [t])[1] ≥ 0 end end + + @testset "Typed tests" begin + ptsP = ManifoldsBase.DefaultMPoint.(pts) + end end end end From c7cdad46aee41ca4ed48973429190e1e3277e3e7 Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Mon, 18 Nov 2019 10:07:02 +0100 Subject: [PATCH 08/12] adds ArrayManifold to Base. --- README.md | 21 ++++- src/ArrayManifold.jl | 172 +++++++++++++++++++++++++++++++++++++++++ src/DefaultManifold.jl | 16 ---- src/ManifoldsBase.jl | 9 ++- 4 files changed, 198 insertions(+), 20 deletions(-) create mode 100644 src/ArrayManifold.jl diff --git a/README.md b/README.md index 11ebfa57..60c01a5a 100644 --- a/README.md +++ b/README.md @@ -4,10 +4,27 @@ Basic interface for manifolds in Julia. +The project [`Manifolds.jl`](https://github.com/JuliaNLSolvers/Manifolds.jl) +is based on this interface and provides a variety of manifolds. + +## Functions on a Manifold + +## DefaultManifold + This interface includes a simple `DefaultManifold`, which is a reduced version of the [`Eucliean`](https://github.com/JuliaNLSolvers/Manifolds.jl/blob/master/src/Euclidean.jl) manifold from [`Manifolds.jl`](https://github.com/JuliaNLSolvers/Manifolds.jl) such that the interface functions can be tested. -The project [`Manifolds.jl`](https://github.com/JuliaNLSolvers/Manifolds.jl) -is based on this interface and provides a variety of manifolds. +## Array Manifold + +The `ArrayManifold` further illustrates, how one can also used types to +represent points on a manifold, tangent vectors as well as cotangent vectors, +where values are encapsulated in a certain type. + +In general this might be used for manifolds, where these three types are represented +by more complex data structures or when it is neccssary to distinguish these +actually by type. + +This adds a semantic layer to the interface and the default implementation of +`ArrayManifold` adds checks to all inputs and outputs of typed data. \ No newline at end of file diff --git a/src/ArrayManifold.jl b/src/ArrayManifold.jl new file mode 100644 index 00000000..8ed674fb --- /dev/null +++ b/src/ArrayManifold.jl @@ -0,0 +1,172 @@ +""" + ArrayManifold{M <: Manifold} <: Manifold + +A manifold to encapsulate manifolds working on array representations of +`MPoints` and `TVectors` in a transparent way, such that for these manifolds its +not necessary to introduce explicit types for the points and tangent vectors, +but they are encapusalted/stripped automatically when needed. +""" +struct ArrayManifold{M <: Manifold} <: Manifold + manifold::M +end +convert(::Type{M},m::ArrayManifold{M}) where M <: Manifold = m.manifold +convert(::Type{ArrayManifold{M}},m::M) where M <: Manifold = ArrayManifold(M) + +manifold_dimension(M::ArrayManifold) = manifold_dimension(M.manifold) + +struct ArrayMPoint{V <: AbstractArray{<:Number}} <: MPoint + value::V +end +convert(::Type{V},x::ArrayMPoint{V}) where V <: AbstractArray{<:Number} = x.value +convert(::Type{ArrayMPoint{V}},x::V) where V <: AbstractArray{<:Number} = ArrayPoint{V}(x) +eltype(::Type{ArrayMPoint{V}}) where V = eltype(V) +similar(x::ArrayMPoint) = ArrayMPoint(similar(x.value)) +similar(x::ArrayMPoint, ::Type{T}) where T = ArrayMPoint(similar(x.value, T)) +function copyto!(x::ArrayMPoint, y::ArrayMPoint) + copyto!(x.value, y.value) + return x +end + +struct ArrayTVector{V <: AbstractArray{<:Number}} <: TVector + value::V +end +convert(::Type{V},v::ArrayTVector{V}) where V <: AbstractArray{<:Number} = v.value +convert(::Type{ArrayTVector{V}},v::V) where V <: AbstractArray{<:Number} = ArrayTVector{V}(v) +eltype(::Type{ArrayTVector{V}}) where V = eltype(V) +similar(x::ArrayTVector) = ArrayTVector(similar(x.value)) +similar(x::ArrayTVector, ::Type{T}) where T = ArrayTVector(similar(x.value, T)) +function copyto!(x::ArrayTVector, y::ArrayTVector) + copyto!(x.value, y.value) + return x +end + +(+)(v1::ArrayTVector, v2::ArrayTVector) = ArrayTVector(v1.value + v2.value) +(-)(v1::ArrayTVector, v2::ArrayTVector) = ArrayTVector(v1.value - v2.value) +(-)(v::ArrayTVector) = ArrayTVector(-v.value) +(*)(a::Number, v::ArrayTVector) = ArrayTVector(a*v.value) + +struct ArrayCoTVector{V <: AbstractArray{<:Number}} <: TVector + value::V +end +convert(::Type{V},v::ArrayCoTVector{V}) where V <: AbstractArray{<:Number} = v.value +convert(::Type{ArrayCoTVector{V}},v::V) where V <: AbstractArray{<:Number} = ArrayCoTVector{V}(v) +eltype(::Type{ArrayCoTVector{V}}) where V = eltype(V) +similar(x::ArrayCoTVector) = ArrayCoTVector(similar(x.value)) +similar(x::ArrayCoTVector, ::Type{T}) where T = ArrayCoTVector(similar(x.value, T)) +function copyto!(x::ArrayCoTVector, y::ArrayCoTVector) + copyto!(x.value, y.value) + return x +end + +(+)(v1::ArrayCoTVector, v2::ArrayCoTVector) = ArrayCoTVector(v1.value + v2.value) +(-)(v1::ArrayCoTVector, v2::ArrayCoTVector) = ArrayCoTVector(v1.value - v2.value) +(-)(v::ArrayCoTVector) = ArrayCoTVector(-v.value) +(*)(a::Number, v::ArrayCoTVector) = ArrayCoTVector(a*v.value) + +array_value(x::AbstractArray) = x +array_value(x::ArrayMPoint) = x.value +array_value(v::ArrayTVector) = v.value +array_value(v::ArrayCoTVector) = v.value + + +function isapprox(M::ArrayManifold, x, y; kwargs...) + is_manifold_point(M, x; kwargs...) + is_manifold_point(M, y; kwargs...) + return isapprox(M.manifold, array_value(x), array_value(y); kwargs...) +end + +function isapprox(M::ArrayManifold, x, v, w; kwargs...) + is_manifold_point(M, x; kwargs...) + is_tangent_vector(M, x, v; kwargs...) + is_tangent_vector(M, x, w; kwargs...) + return isapprox(M.manifold, array_value(x), array_value(v), array_value(w); kwargs...) +end + +function project_tangent!(M::ArrayManifold, w, x, v; kwargs...) + is_manifold_point(M, x; kwargs...) + project_tangent!(M.manifold, w.value, array_value(x), array_value(v)) + is_tangent_vector(M, x, w; kwargs...) + return w +end + +function distance(M::ArrayManifold, x, y; kwargs...) + is_manifold_point(M, x; kwargs...) + is_manifold_point(M, y; kwargs...) + return distance(M.manifold, array_value(x), array_value(y)) +end + +function inner(M::ArrayManifold, x, v, w; kwargs...) + is_manifold_point(M, x; kwargs...) + is_tangent_vector(M, x, v; kwargs...) + is_tangent_vector(M, x, w; kwargs...) + return inner(M.manifold, array_value(x), array_value(v), array_value(w)) +end + +function exp(M::ArrayManifold, x, v; kwargs...) + is_manifold_point(M, x; kwargs...) + is_tangent_vector(M, x, v; kwargs...) + y = ArrayMPoint(exp(M.manifold, array_value(x), array_value(v))) + is_manifold_point(M, y; kwargs...) + return y +end + +function exp!(M::ArrayManifold, y, x, v; kwargs...) + is_manifold_point(M, x; kwargs...) + is_tangent_vector(M, x, v; kwargs...) + exp!(M.manifold, array_value(y), array_value(x), array_value(v)) + is_manifold_point(M, y; kwargs...) + return y +end + +function log(M::ArrayManifold, x, y; kwargs...) + is_manifold_point(M, x; kwargs...) + is_manifold_point(M, y; kwargs...) + v = ArrayTVector(log(M.manifold, array_value(x), array_value(y))) + is_tangent_vector(M, x, v; kwargs...) + return v +end + +function log!(M::ArrayManifold, v, x, y; kwargs...) + is_manifold_point(M, x; kwargs...) + is_manifold_point(M, y; kwargs...) + log!(M.manifold, array_value(v), array_value(x), array_value(y)) + is_tangent_vector(M, x, v; kwargs...) + return v +end + +function zero_tangent_vector!(M::ArrayManifold, v, x; kwargs...) + is_manifold_point(M, x; kwargs...) + zero_tangent_vector!(M.manifold, array_value(v), array_value(x); kwargs...) + is_tangent_vector(M, x, v; kwargs...) + return v +end + +function zero_tangent_vector(M::ArrayManifold, x; kwargs...) + is_manifold_point(M, x; kwargs...) + w = zero_tangent_vector(M.manifold, array_value(x)) + is_tangent_vector(M, x, w; kwargs...) + return w +end + +function vector_transport_to(M::ArrayManifold, x, v, y) + return vector_transport_to(M.manifold, + array_value(x), + array_value(v), + array_value(y)) +end + +function vector_transport_to!(M::ArrayManifold, vto, x, v, y) + return vector_transport_to!(M.manifold, + array_value(vto), + array_value(x), + array_value(v), + array_value(y)) +end + +function is_manifold_point(M::ArrayManifold, x::MPoint; kwargs...) + return is_manifold_point(M.manifold, array_value(x); kwargs...) +end + +function is_tangent_vector(M::ArrayManifold, x::MPoint, v::TVector; kwargs...) + return is_tangent_vector(M.manifold, array_value(x), array_value(v); kwargs...) +end diff --git a/src/DefaultManifold.jl b/src/DefaultManifold.jl index c5c681b2..102b0910 100644 --- a/src/DefaultManifold.jl +++ b/src/DefaultManifold.jl @@ -13,22 +13,6 @@ involved variabes. struct DefaultManifold{T<:Tuple} <: Manifold where {T} end DefaultManifold(n::Vararg{Int,N}) where N = DefaultManifold{Tuple{n...}}() -struct DefaultMPoint{T} <: MPoint where {T} - data::T -end -struct DefaultTVector{T} <: TVector where {T} - data::T -end -struct DefaultCoTVector{T} <: CoTVector where {T} - data::T -end -for T in [DefaultMPoint,DefaultTVector,DefaultCoTVector] - Base.eltype(p::T) = eltype(p.data) - Base.similar(p::T, ::Type{EltypeNew}) where EltypeNew = T(similar(p.data, EltypeNew)) - Base.similar(p::T) = T(similar(p.data)) - Base.copy(p::T) = T(copy(p.data)) -end - @generated representation_size(::DefaultManifold{T}) where {T} = Tuple(T.parameters...) @generated manifold_dimension(::DefaultManifold{T}) where {T} = *(T.parameters...) @inline inner(::DefaultManifold, x, v, w) = dot(v, w) diff --git a/src/ManifoldsBase.jl b/src/ManifoldsBase.jl index 3c5360f3..4f65e86e 100644 --- a/src/ManifoldsBase.jl +++ b/src/ManifoldsBase.jl @@ -644,14 +644,19 @@ function check_tangent_vector(M::Manifold, x, v; kwargs...) end end - +include("ArrayManifold.jl") include("DefaultManifold.jl") export Manifold, MPoint, TVector, CoTVector, - ParallelTransport, + ArrayManifold, + ArrayMPoint, + ArrayTVector, + ArrayCoTVector + +export ParallelTransport, ProjectionTransport export base_manifold, From 9db1829ddb393ff717e44f0008cc51330cf249db Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Mon, 18 Nov 2019 10:16:56 +0100 Subject: [PATCH 09/12] Adds documentation to the `ArrayManifold` types and new functions as well as testing. --- src/ArrayManifold.jl | 30 +++++++++++++++++++++++++++++- test/array_manifold.jl | 26 ++++++++++++++++++++++++++ test/default_manifold.jl | 4 ---- test/runtests.jl | 1 + 4 files changed, 56 insertions(+), 5 deletions(-) create mode 100644 test/array_manifold.jl diff --git a/src/ArrayManifold.jl b/src/ArrayManifold.jl index 8ed674fb..bf30f7c6 100644 --- a/src/ArrayManifold.jl +++ b/src/ArrayManifold.jl @@ -4,7 +4,7 @@ A manifold to encapsulate manifolds working on array representations of `MPoints` and `TVectors` in a transparent way, such that for these manifolds its not necessary to introduce explicit types for the points and tangent vectors, -but they are encapusalted/stripped automatically when needed. +but they are encapsulated/stripped automatically when needed. """ struct ArrayManifold{M <: Manifold} <: Manifold manifold::M @@ -14,6 +14,13 @@ convert(::Type{ArrayManifold{M}},m::M) where M <: Manifold = ArrayManifold(M) manifold_dimension(M::ArrayManifold) = manifold_dimension(M.manifold) +""" + ArrayMPoint <: MPoint + +represent a point on an [`ArrayManfold`](@ref), i.e. on a manifold where data +can be represented by arrays. The array is stored internally and semantically +this distinguished the value from [`ArrayTVector`](@ref)s and [`ArrayCoTVector`](@ref)s +""" struct ArrayMPoint{V <: AbstractArray{<:Number}} <: MPoint value::V end @@ -27,6 +34,13 @@ function copyto!(x::ArrayMPoint, y::ArrayMPoint) return x end +""" + ArrayTVector <: TVector + +represent a tangent vector an [`ArrayManfold`](@ref), i.e. on a manifold where data +can be represented by arrays. The array is stored internally and semantically +this distinguished the value from [`ArrayMPoint`](@ref)s and [`ArrayCoTVector`](@ref)s +""" struct ArrayTVector{V <: AbstractArray{<:Number}} <: TVector value::V end @@ -45,6 +59,13 @@ end (-)(v::ArrayTVector) = ArrayTVector(-v.value) (*)(a::Number, v::ArrayTVector) = ArrayTVector(a*v.value) +""" + ArrayCoTVector <: CoTVector + +represent a cotangent vector an [`ArrayManfold`](@ref), i.e. on a manifold where data +can be represented by arrays. The array is stored internally and semantically +this distinguished the value from [`ArrayMPoint`](@ref)s and [`ArrayTVector`](@ref)s +""" struct ArrayCoTVector{V <: AbstractArray{<:Number}} <: TVector value::V end @@ -63,6 +84,13 @@ end (-)(v::ArrayCoTVector) = ArrayCoTVector(-v.value) (*)(a::Number, v::ArrayCoTVector) = ArrayCoTVector(a*v.value) +""" + array_value(x) + +returns the internal array value of a [`ArrayMPoint`](@ref), [`ArrayTVector`](@ref) +or [`ArrayCoTVector`](@ref) if the value `x` is encapsulated as such, otherwise +if `x` is already an array, it just returns `x` +""" array_value(x::AbstractArray) = x array_value(x::ArrayMPoint) = x.value array_value(v::ArrayTVector) = v.value diff --git a/test/array_manifold.jl b/test/array_manifold.jl new file mode 100644 index 00000000..597971fe --- /dev/null +++ b/test/array_manifold.jl @@ -0,0 +1,26 @@ +using ManifoldsBase +using LinearAlgebra + +@testset "Array manifold" begin + M = ManifoldsBase.DefaultManifold(3) + A = ArrayManifold(M) + x = [1., 0., 0.] + y = 1/sqrt(2)*[1., 1., 0.] + z = [0., 1., 0.] + v = log(M,x,y) + x2 = ArrayMPoint(x) + y2 = ArrayMPoint(y) + v2 = log(A,x,y) # auto convert + y2 = exp(A,x,v2) + w = log(M,x,z) + w2 = log(A,x,z; atol=10^(-15)) + @test isapprox(y2.value,y) + @test distance(A,x,y) == distance(M,x,y) + @test norm(A,x,v) == norm(M,x,v) + @test inner(A,x,v2,w2; atol=10^(-15)) == inner(M,x,v,w) + @test isapprox(A, x2, y2 ) == isapprox(M,x,y) + @test isapprox(A,x,y) == isapprox(A,x2,y2) + @test isapprox(A,x, v2,v2 ) == isapprox(M,x,v,v) + @test isapprox(A, exp(A,x,v),y2) + @test isapprox(A, zero_tangent_vector(A,x), zero_tangent_vector(M,x)) +end diff --git a/test/default_manifold.jl b/test/default_manifold.jl index 6a4a56ae..ed7f01b6 100644 --- a/test/default_manifold.jl +++ b/test/default_manifold.jl @@ -158,10 +158,6 @@ using Test @test ReverseDiff.gradient(retract_f, [t])[1] ≥ 0 end end - - @testset "Typed tests" begin - ptsP = ManifoldsBase.DefaultMPoint.(pts) - end end end end diff --git a/test/runtests.jl b/test/runtests.jl index 33a646cc..983041a4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,3 +1,4 @@ include("empty_manifold.jl") include("default_manifold.jl") +include("array_manifold.jl") From ec5005216a82a793c0cddbefbfcb2a533bae744b Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Mon, 18 Nov 2019 10:37:50 +0100 Subject: [PATCH 10/12] adapt vector_transport to the new interface in ArrayManifold. --- src/ArrayManifold.jl | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/src/ArrayManifold.jl b/src/ArrayManifold.jl index bf30f7c6..9aeaa417 100644 --- a/src/ArrayManifold.jl +++ b/src/ArrayManifold.jl @@ -176,19 +176,22 @@ function zero_tangent_vector(M::ArrayManifold, x; kwargs...) return w end -function vector_transport_to(M::ArrayManifold, x, v, y) - return vector_transport_to(M.manifold, - array_value(x), - array_value(v), - array_value(y)) +function vector_transport_to!(M::ArrayManifold, vto, x, v, y, m::AbstractVectorTransportMethod) + return vector_transport_to!(M.manifold, + array_value(vto), + array_value(x), + array_value(v), + array_value(y), + m) end -function vector_transport_to!(M::ArrayManifold, vto, x, v, y) - return vector_transport_to!(M.manifold, +function vector_transport_along!(M::ArrayManifold, vto, x, v, c, m::AbstractVectorTransportMethod) + return vector_transport_along!(M.manifold, array_value(vto), array_value(x), array_value(v), - array_value(y)) + c, + m) end function is_manifold_point(M::ArrayManifold, x::MPoint; kwargs...) From 59c8efa8a2f4a2f370ecbbed629c29c313cdc130 Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Mon, 18 Nov 2019 15:24:08 +0100 Subject: [PATCH 11/12] fixes a few minor bugs and increases code coverage. --- src/ArrayManifold.jl | 4 +-- src/ManifoldsBase.jl | 4 +++ test/array_manifold.jl | 62 ++++++++++++++++++++++++++++++++++++------ 3 files changed, 59 insertions(+), 11 deletions(-) diff --git a/src/ArrayManifold.jl b/src/ArrayManifold.jl index 9aeaa417..bc85e219 100644 --- a/src/ArrayManifold.jl +++ b/src/ArrayManifold.jl @@ -10,7 +10,7 @@ struct ArrayManifold{M <: Manifold} <: Manifold manifold::M end convert(::Type{M},m::ArrayManifold{M}) where M <: Manifold = m.manifold -convert(::Type{ArrayManifold{M}},m::M) where M <: Manifold = ArrayManifold(M) +convert(::Type{ArrayManifold{M}},m::M) where M <: Manifold = ArrayManifold(m) manifold_dimension(M::ArrayManifold) = manifold_dimension(M.manifold) @@ -25,7 +25,7 @@ struct ArrayMPoint{V <: AbstractArray{<:Number}} <: MPoint value::V end convert(::Type{V},x::ArrayMPoint{V}) where V <: AbstractArray{<:Number} = x.value -convert(::Type{ArrayMPoint{V}},x::V) where V <: AbstractArray{<:Number} = ArrayPoint{V}(x) +convert(::Type{ArrayMPoint{V}},x::V) where V <: AbstractArray{<:Number} = ArrayMPoint{V}(x) eltype(::Type{ArrayMPoint{V}}) where V = eltype(V) similar(x::ArrayMPoint) = ArrayMPoint(similar(x.value)) similar(x::ArrayMPoint, ::Type{T}) where T = ArrayMPoint(similar(x.value, T)) diff --git a/src/ManifoldsBase.jl b/src/ManifoldsBase.jl index 4f65e86e..53c95da9 100644 --- a/src/ManifoldsBase.jl +++ b/src/ManifoldsBase.jl @@ -3,6 +3,8 @@ module ManifoldsBase import Base: isapprox, exp, log, + convert, + copyto!, angle, eltype, similar, @@ -686,6 +688,8 @@ export base_manifold, retract, retract!, tangent_vector_error, + vector_transport_along, + vector_transport_along!, vector_transport_direction, vector_transport_direction!, vector_transport_to, diff --git a/test/array_manifold.jl b/test/array_manifold.jl index 597971fe..e66f3481 100644 --- a/test/array_manifold.jl +++ b/test/array_manifold.jl @@ -14,13 +14,57 @@ using LinearAlgebra y2 = exp(A,x,v2) w = log(M,x,z) w2 = log(A,x,z; atol=10^(-15)) - @test isapprox(y2.value,y) - @test distance(A,x,y) == distance(M,x,y) - @test norm(A,x,v) == norm(M,x,v) - @test inner(A,x,v2,w2; atol=10^(-15)) == inner(M,x,v,w) - @test isapprox(A, x2, y2 ) == isapprox(M,x,y) - @test isapprox(A,x,y) == isapprox(A,x2,y2) - @test isapprox(A,x, v2,v2 ) == isapprox(M,x,v,v) - @test isapprox(A, exp(A,x,v),y2) - @test isapprox(A, zero_tangent_vector(A,x), zero_tangent_vector(M,x)) + @testset "Types and Conversion" begin + @test convert(typeof(M), A) == M + @test convert(typeof(A),M) == A + + for T in [ArrayMPoint, ArrayTVector, ArrayCoTVector] + p = T(x) + @test convert(typeof(x),p) == x + @test convert(typeof(p),y) == T(y) + @test eltype(p) == eltype(x) + @test typeof(similar(p)) == typeof(p) + @test typeof(similar(p,eltype(x))) == typeof(p) + q = similar(p) + copyto!(q,p) + @test isapprox(A,q,p) + @test ManifoldsBase.array_value(p) == x + @test ManifoldsBase.array_value(x) == x + end + end + @testset "Vector functions" begin + for T in [ArrayTVector, ArrayCoTVector] + a = T(v) + b = T(w) + @test isapprox(A, a+b, T(v+w)) + @test isapprox(A, (a-b), T(v-w) ) + @test isapprox(A, -b, T(-w) ) + @test isapprox(A, 2*a, T(2 .* v) ) + end + end + @testset "Manifold functions" begin + @test manifold_dimension(A) == manifold_dimension(M) + @test isapprox(y2.value,y) + @test distance(A,x,y) == distance(M,x,y) + @test norm(A,x,v) == norm(M,x,v) + @test inner(A,x,v2,w2; atol=10^(-15)) == inner(M,x,v,w) + @test isapprox(A, x2, y2 ) == isapprox(M,x,y) + @test isapprox(A,x,y) == isapprox(A,x2,y2) + @test isapprox(A,x, v2,v2 ) == isapprox(M,x,v,v) + v2s = similar(v2) + project_tangent!(A,v2s,x2,v2) + @test isapprox(A, v2, v2s) + y2s = similar(y2) + exp!(A,y2s,x2,v2) + @test isapprox(A,y2s,y2) + log!(A, v2s, x, y) + @test isapprox(A, x, v2s, v2) + @test isapprox(A, exp(A,x,v),y2) + @test isapprox(A, zero_tangent_vector(A,x), zero_tangent_vector(M,x)) + vector_transport_to!(A, v2s, x2, v2, y2) + @test isapprox(A, x2, v2, v2s) + zero_tangent_vector!(A, v2s, x) + @test isapprox(A, v2s, zero_tangent_vector(M,x)) + @test_throws ErrorException vector_transport_along!(A,v2s,x2,v2,ParallelTransport()) + end end From 611a04b0edcd0c4aad498c843c2386e79dcfc4e8 Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Mon, 18 Nov 2019 15:47:13 +0100 Subject: [PATCH 12/12] adds another eltype test for ArrayManifold point/vector types. --- test/array_manifold.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/array_manifold.jl b/test/array_manifold.jl index e66f3481..08b3069b 100644 --- a/test/array_manifold.jl +++ b/test/array_manifold.jl @@ -22,6 +22,7 @@ using LinearAlgebra p = T(x) @test convert(typeof(x),p) == x @test convert(typeof(p),y) == T(y) + @test eltype(typeof(p)) == eltype(x) @test eltype(p) == eltype(x) @test typeof(similar(p)) == typeof(p) @test typeof(similar(p,eltype(x))) == typeof(p)