diff --git a/Project.toml b/Project.toml index 1401bc7021..8adb2ee620 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.5.8" +version = "0.5.9" [deps] Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" diff --git a/docs/make.jl b/docs/make.jl index 70605f93d1..575a8f5c00 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -64,8 +64,9 @@ makedocs( "Vector bundle" => "manifolds/vector_bundle.md", ], "Manifold decorators" => [ - "Metric manifold" => "manifolds/metric.md", + "Connection manifold" => "manifolds/connection.md", "Group manifold" => "manifolds/group.md", + "Metric manifold" => "manifolds/metric.md", ], ], "Features on Manifolds" => [ diff --git a/docs/src/manifolds/connection.md b/docs/src/manifolds/connection.md new file mode 100644 index 0000000000..798e06395b --- /dev/null +++ b/docs/src/manifolds/connection.md @@ -0,0 +1,39 @@ +# Connection manifold + +A connection manifold always consists of a [topological manifold](https://en.wikipedia.org/wiki/Topological_manifold) together with a [connection](https://en.wikipedia.org/wiki/Connection_(mathematics)) $\Gamma$. + +However, often there is an implicitly assumed (default) connection, like the [`LeviCivitaConnection`](@ref) connection on a Riemannian manifold. +It is not necessary to use this decorator if you implement just one (or the first) connection. +If you later introduce a second, the old (first) connection can be used with the (non [`AbstractConnectionManifold`](@ref)) [`AbstractManifold`](@ref), i.e. without an explicitly stated connection. + +This manifold decorator serves two purposes: + +1. to implement different connections (e.g. in closed form) for one [`AbstractManifold`](@ref) +2. to provide a way to compute geodesics on manifolds, where this [`AbstractAffineConnection`](@ref) does not yield a closed formula. + +An example of usage can be found in Cartan-Schouten connections, see [`AbstractCartanSchoutenConnection`](@ref). + +```@contents +Pages = ["connection.md"] +Depth = 2 +``` + +## Types + +```@autodocs +Modules = [Manifolds, ManifoldsBase] +Pages = ["manifolds/ConnectionManifold.jl"] +Order = [:type] +``` + +## Functions + +```@autodocs +Modules = [Manifolds, ManifoldsBase] +Pages = ["manifolds/ConnectionManifold.jl"] +Order = [:function] +``` + +## [Charts and bases of vector spaces](@id connections_charts) + +All connection-related functions take a basis of a vector space as one of the arguments. This is needed because generally there is no way to define these functions without referencing a basis. In some cases there is no need to be explicit about this basis, and then for example a [`DefaultOrthonormalBasis`](@ref) object can be used. In cases where being explicit about these bases is needed, for example when using multiple charts, a basis can be specified, for example using [`induced_basis`](@ref Main.Manifolds.induced_basis). diff --git a/docs/src/manifolds/group.md b/docs/src/manifolds/group.md index 12b1d280fe..5420d7903e 100644 --- a/docs/src/manifolds/group.md +++ b/docs/src/manifolds/group.md @@ -4,7 +4,10 @@ Lie groups, groups that are [`AbstractManifold`](@ref)s with a smooth binary gro The common addition and multiplication group operations of [`AdditionOperation`](@ref) and [`MultiplicationOperation`](@ref) are provided, though their behavior may be customized for a specific group. +There are short introductions at the beginning of each subsection. They briefly mention what is available with links to more detailed descriptions. + #### Contents + ```@contents Pages = ["group.md"] Depth = 3 @@ -12,8 +15,17 @@ Depth = 3 ## Groups +The following operations are available for group manifolds: + +* [`identity`](@ref): get the identity of the group. +* [`inv`](@ref): get the inverse of a given element. +* [`compose`](@ref): compose two given elements of a group. + ### Group manifold +[`GroupManifold`](@ref) adds a group structure to the wrapped manifold. +It does not affect metric (or connection) structure of the wrapped manifold, however it can to be further wrapped in [`MetricManifold`](@ref) to get invariant metrics, or in a [`ConnectionManifold`](@ref) to equip it with a Cartan-Schouten connection. + ```@autodocs Modules = [Manifolds] Pages = ["groups/group.jl"] @@ -86,6 +98,29 @@ Order = [:type, :function] ## Group actions +Group actions represent actions of a given group on a specified manifold. +The following operations are available: + +* [`apply`](@ref): performs given action of an element of the group on an object of compatible type. +* [`apply_diff`](@ref): differential of [`apply`](@ref) with respect to the object it acts upon. +* [`direction`](@ref): tells whether a given action is [`LeftAction`](@ref) or [`RightAction`](@ref). +* [`inverse_apply`](@ref): performs given action of the inverse of an element of the group on an object of compatible type. By default inverts the element and calls [`apply`](@ref) but it may be have a faster implementation for some actions. +* [`inverse_apply_diff`](@ref): counterpart of [`apply_diff`](@ref) for [`inverse_apply`](@ref). +* [`optimal_alignment`](@ref): determine the element of a group that, when it acts upon a point, produces the element closest to another given point in the metric of the G-manifold. + +Furthermore, group operation action features the following: + +* [`translate`](@ref Main.Manifolds.translate): an operation that performs either left ([`LeftAction`](@ref)) or right ([`RightAction`](@ref)) translation. This is by default performed by calling [`compose`](@ref) with appropriate order of arguments. This function is separated from `compose` mostly to easily represent its differential, [`translate_diff`](@ref). +* [`translate_diff`](@ref): differential of [`translate`](@ref Main.Manifolds.translate) with respect to the point being translated. +* [`adjoint_action`](@ref): adjoint action of a given element of a Lie group on an element of its Lie algebra. +* [`lie_bracket`](@ref): Lie bracket of two vectors from a Lie algebra corresponding to a given group. + +The following group actions are available: + +* Group operation action [`GroupOperationAction`](@ref) that describes action of a group on itself. +* [`RotationAction`](@ref), that is action of [`SpecialOrthogonal`](@ref) group on different manifolds. +* [`TranslationAction`](@ref), which is the action of [`TranslationGroup`](@ref) group on different manifolds. + ```@autodocs Modules = [Manifolds] Pages = ["groups/group_action.jl"] @@ -123,3 +158,11 @@ Modules = [Manifolds] Pages = ["groups/metric.jl"] Order = [:type, :function] ``` + +## Cartan-Schouten connections + +```@autodocs +Modules = [Manifolds] +Pages = ["groups/connections.jl"] +Order = [:type, :function] +``` diff --git a/docs/src/manifolds/metric.md b/docs/src/manifolds/metric.md index 626ab04a45..dd34512be0 100644 --- a/docs/src/manifolds/metric.md +++ b/docs/src/manifolds/metric.md @@ -17,6 +17,8 @@ Pages = ["metric.md"] Depth = 2 ``` +Note that a metric manifold is an [`AbstractConnectionManifold`](@ref) with the [`LeviCivitaConnection`](@ref) of the metric $g$, and thus a large part of metric manifold's functionality relies on this. + Let's first look at the provided types. ## Types @@ -47,6 +49,6 @@ Order = [:function] ## Metrics, charts and bases of vector spaces -All metric-related functions take a basis of a vector space as one of the arguments. This needed because generally there is no way to define these functions without referencing a basis. In some cases there is no need to be explicit about this basis, and then for example a [`DefaultOrthonormalBasis`](@ref) object can be used. In cases where being explicit about these bases is needed, for example when using multiple charts, a basis can be specified, for example using [`induced_basis`](@ref Main.Manifolds.induced_basis). +Metric-related functions, similarly to connection-related functions, need to operate in a basis of a vector space, see [here](@ref connections_charts). Metric-related functions can take bases of associated tangent spaces as arguments. For example [`local_metric`](@ref) can take the basis of the tangent space it is supposed to operate on instead of a custom basis of the space of symmetric bilinear operators. diff --git a/docs/src/misc/notation.md b/docs/src/misc/notation.md index 9c05384b14..88e1a810d0 100644 --- a/docs/src/misc/notation.md +++ b/docs/src/misc/notation.md @@ -10,6 +10,7 @@ Within the documented functions, the utf8 symbols are used whenever possible, as | Symbol | Description | Also used | Comment | |:--:|:--------------- |:--:|:-- | | ``\tau_p`` | action map by group element ``p`` | ``\mathrm{L}_p``, ``\mathrm{R}_p`` | either left or right | +| ``\operatorname{Ad}_p(X)`` | adjoint action of element ``p`` of a Lie group on the element ``X`` of the corresponding Lie algebra | | | | ``\times`` | Cartesian product of two manifolds | | see [`ProductManifold`](@ref) | | ``^{\wedge}`` | (n-ary) Cartesian power of a manifold | | see [`PowerManifold`](@ref) | | ``a`` | coordinates of a point in a chart | | see [`get_parameters`](@ref) | diff --git a/src/Manifolds.jl b/src/Manifolds.jl index a62acbaf7a..8e7e5826e5 100644 --- a/src/Manifolds.jl +++ b/src/Manifolds.jl @@ -149,6 +149,8 @@ include("manifolds/VectorBundle.jl") include("distributions.jl") include("projected_distribution.jl") +include("manifolds/ConnectionManifold.jl") + # It's included early to ensure visibility of `Identity` include("groups/group.jl") @@ -195,6 +197,7 @@ include("manifolds/Multinomial.jl") include("manifolds/Oblique.jl") include("manifolds/EssentialManifold.jl") +include("groups/connections.jl") include("groups/metric.jl") include("groups/group_action.jl") include("groups/group_operation_action.jl") @@ -352,6 +355,10 @@ export AbstractVectorTransportMethod, export PoleLadderTransport, SchildsLadderTransport export PowerVectorTransport, ProductVectorTransport export AbstractEmbeddedManifold +export AbstractAffineConnection, + AbstractConnectionManifold, ConnectionManifold, LeviCivitaConnection +export AbstractCartanSchoutenConnection, + CartanSchoutenMinus, CartanSchoutenPlus, CartanSchoutenZero export AbstractMetric, RiemannianMetric, LorentzMetric, @@ -528,7 +535,9 @@ export AbstractGroupAction, SpecialOrthogonal, TranslationGroup, TranslationAction -export affine_matrix, +export adjoint_action, + adjoint_action!, + affine_matrix, apply, apply!, apply_diff, @@ -559,6 +568,8 @@ export affine_matrix, inverse_translate!, inverse_translate_diff, inverse_translate_diff!, + lie_bracket, + lie_bracket!, make_identity, optimal_alignment, optimal_alignment!, diff --git a/src/groups/array_manifold.jl b/src/groups/array_manifold.jl index 6c91b5da9f..2f9c2411b4 100644 --- a/src/groups/array_manifold.jl +++ b/src/groups/array_manifold.jl @@ -4,6 +4,24 @@ array_point(p) = ValidationMPoint(p) array_point(p::ValidationMPoint) = p array_point(e::Identity) = Identity(e.group, array_point(e.p)) +function adjoint_action(M::ValidationManifold, p, X; kwargs...) + is_point(M, p, true; kwargs...) + eM = make_identity(M.manifold, array_value(p)) + is_vector(M, eM, X, true; kwargs...) + Y = ValidationTVector(adjoint_action(M.manifold, array_value(p), array_value(X))) + is_vector(M, eM, Y, true; kwargs...) + return Y +end + +function adjoint_action!(M::ValidationManifold, Y, p, X; kwargs...) + is_point(M, p, true; kwargs...) + eM = make_identity(M.manifold, array_value(p)) + is_vector(M, eM, X, true; kwargs...) + adjoint_action!(M.manifold, array_value(Y), array_value(p), array_value(X)) + is_vector(M, eM, Y, true; kwargs...) + return Y +end + function Base.inv(M::ValidationManifold, p; kwargs...) is_point(M, p, true; kwargs...) q = array_point(inv(M.manifold, array_value(p))) @@ -32,6 +50,24 @@ function identity!(M::ValidationManifold, q, p; kwargs...) return q end +function lie_bracket(M::ValidationManifold, X, Y) + eM = make_identity(M.manifold, array_value(X)) + is_vector(M, eM, X, true) + is_vector(M, eM, Y, true) + Z = ValidationTVector(lie_bracket(M.manifold, array_value(X), array_value(Y))) + is_vector(M, eM, Z, true) + return Z +end + +function lie_bracket!(M::ValidationManifold, Z, X, Y) + eM = make_identity(M.manifold, array_value(X)) + is_vector(M, eM, X, true) + is_vector(M, eM, Y, true) + lie_bracket!(M.manifold, array_value(Z), array_value(X), array_value(Y)) + is_vector(M, eM, Z, true) + return Z +end + function compose(M::ValidationManifold, p, q; kwargs...) is_point(M, p, true; kwargs...) is_point(M, q, true; kwargs...) diff --git a/src/groups/circle_group.jl b/src/groups/circle_group.jl index 46e46265f0..2f82394086 100644 --- a/src/groups/circle_group.jl +++ b/src/groups/circle_group.jl @@ -14,6 +14,10 @@ invariant_metric_dispatch(::CircleGroup, ::ActionDirection) = Val(true) default_metric_dispatch(::MetricManifold{ℂ,CircleGroup,EuclideanMetric}) = Val(true) +adjoint_action(::CircleGroup, p, X) = X + +adjoint_action!(::CircleGroup, Y, p, X) = copyto!(Y, X) + function compose(G::CircleGroup, p::AbstractVector, q::AbstractVector) return map(compose, repeated(G), p, q) end @@ -46,6 +50,10 @@ function inverse_translate( return map(/, q, p) end +lie_bracket(::CircleGroup, X, Y) = zero(X) + +lie_bracket!(::CircleGroup, Z, X, Y) = fill!(Z, 0) + translate_diff(::GT, p, q, X, ::ActionDirection) where {GT<:CircleGroup} = map(*, p, X) function translate_diff( ::GT, diff --git a/src/groups/connections.jl b/src/groups/connections.jl new file mode 100644 index 0000000000..8fcee83d33 --- /dev/null +++ b/src/groups/connections.jl @@ -0,0 +1,182 @@ + +""" + AbstractCartanSchoutenConnection + +Abstract type for Cartan-Schouten connections, that is connections whose geodesics +going through group identity are one-parameter subgroups. See [^Pennec2020] for details. + +[^Pennec2020]: + > X. Pennec and M. Lorenzi, “5 - Beyond Riemannian geometry: The affine connection + > setting for transformation groups,” in Riemannian Geometric Statistics in Medical Image + > Analysis, X. Pennec, S. Sommer, and T. Fletcher, Eds. Academic Press, 2020, pp. 169–229. + > doi: 10.1016/B978-0-12-814725-2.00012-1. +""" +abstract type AbstractCartanSchoutenConnection <: AbstractAffineConnection end + +""" + CartanSchoutenMinus + +The unique Cartan-Schouten connection such that all left-invariant vector fields are +globally defined by their value at identity. It is biinvariant with respect to the group +operation. +""" +struct CartanSchoutenMinus <: AbstractCartanSchoutenConnection end + +""" + CartanSchoutenPlus + +The unique Cartan-Schouten connection such that all right-invariant vector fields are +globally defined by their value at identity. It is biinvariant with respect to the group +operation. +""" +struct CartanSchoutenPlus <: AbstractCartanSchoutenConnection end + +""" + CartanSchoutenZero + +The unique torsion-free Cartan-Schouten connection. It is biinvariant with respect to the +group operation. + +If the metric on the underlying manifold is bi-invariant then it is equivalent to the +Levi-Civita connection of that metric. +""" +struct CartanSchoutenZero <: AbstractCartanSchoutenConnection end + +const CartanSchoutenMinusGroup{𝔽,M} = ConnectionManifold{𝔽,M,CartanSchoutenMinus} +const CartanSchoutenPlusGroup{𝔽,M} = ConnectionManifold{𝔽,M,CartanSchoutenPlus} +const CartanSchoutenZeroGroup{𝔽,M} = ConnectionManifold{𝔽,M,CartanSchoutenZero} + +""" + exp!(M::ConnectionManifold{𝔽,<:AbstractGroupManifold{𝔽},<:AbstractCartanSchoutenConnection}, q, p, X) where {𝔽} + +Compute the exponential map on the [`ConnectionManifold`](@ref) `M` with a Cartan-Schouten +connection. See Sections 5.3.2 and 5.3.3 of [^Pennec2020] for details. + +[^Pennec2020]: + > X. Pennec and M. Lorenzi, “5 - Beyond Riemannian geometry: The affine connection + > setting for transformation groups,” in Riemannian Geometric Statistics in Medical Image + > Analysis, X. Pennec, S. Sommer, and T. Fletcher, Eds. Academic Press, 2020, pp. 169–229. + > doi: 10.1016/B978-0-12-814725-2.00012-1. +""" +function exp!( + M::ConnectionManifold{𝔽,<:AbstractGroupManifold{𝔽},<:AbstractCartanSchoutenConnection}, + q, + p, + X, +) where {𝔽} + Y = inverse_translate_diff(M.manifold, q, p, X) + return compose!(M.manifold, q, p, group_exp(M.manifold, Y)) +end + +""" + log!(M::ConnectionManifold{𝔽,<:AbstractGroupManifold{𝔽},<:AbstractCartanSchoutenConnection}, Y, p, q) where {𝔽} + +Compute the logarithmic map on the [`ConnectionManifold`](@ref) `M` with a Cartan-Schouten +connection. See Sections 5.3.2 and 5.3.3 of [^Pennec2020] for details. + +[^Pennec2020]: + > X. Pennec and M. Lorenzi, “5 - Beyond Riemannian geometry: The affine connection + > setting for transformation groups,” in Riemannian Geometric Statistics in Medical Image + > Analysis, X. Pennec, S. Sommer, and T. Fletcher, Eds. Academic Press, 2020, pp. 169–229. + > doi: 10.1016/B978-0-12-814725-2.00012-1. +""" +function log!( + M::ConnectionManifold{𝔽,<:AbstractGroupManifold{𝔽},<:AbstractCartanSchoutenConnection}, + Y, + p, + q, +) where {𝔽} + pinvq = compose(M.manifold, inv(M.manifold, p), q) + group_log!(M.manifold, Y, pinvq) + return translate_diff!(M.manifold, Y, p, Identity(M.manifold, p), Y) +end + +""" + vector_transport_to(M::CartanSchoutenMinusGroup, p, X, q, ::ParallelTransport) + +Transport tangent vector `X` at point `p` on the group manifold `M` with the +[`CartanSchoutenMinus`](@ref) connection to point `q`. See [^Pennec2020] for details. + +[^Pennec2020]: + > X. Pennec and M. Lorenzi, “5 - Beyond Riemannian geometry: The affine connection + > setting for transformation groups,” in Riemannian Geometric Statistics in Medical Image + > Analysis, X. Pennec, S. Sommer, and T. Fletcher, Eds. Academic Press, 2020, pp. 169–229. + > doi: 10.1016/B978-0-12-814725-2.00012-1. +""" +vector_transport_to(M::CartanSchoutenMinusGroup, p, X, q, ::ParallelTransport) + +function vector_transport_to!(M::CartanSchoutenMinusGroup, Y, p, X, q, ::ParallelTransport) + return inverse_translate_diff!(M.manifold, Y, q, p, X, LeftAction()) +end + +""" + vector_transport_to(M::CartanSchoutenPlusGroup, p, X, q, ::ParallelTransport) + +Transport tangent vector `X` at point `p` on the group manifold `M` with the +[`CartanSchoutenPlus`](@ref) connection to point `q`. See [^Pennec2020] for details. + +[^Pennec2020]: + > X. Pennec and M. Lorenzi, “5 - Beyond Riemannian geometry: The affine connection + > setting for transformation groups,” in Riemannian Geometric Statistics in Medical Image + > Analysis, X. Pennec, S. Sommer, and T. Fletcher, Eds. Academic Press, 2020, pp. 169–229. + > doi: 10.1016/B978-0-12-814725-2.00012-1. +""" +vector_transport_to(M::CartanSchoutenPlusGroup, p, X, q, ::ParallelTransport) + +function vector_transport_to!(M::CartanSchoutenPlusGroup, Y, p, X, q, ::ParallelTransport) + return inverse_translate_diff!(M.manifold, Y, q, p, X, RightAction()) +end + +""" + vector_transport_direction(M::CartanSchoutenZeroGroup, ::Identity, X, d, ::ParallelTransport) + +Transport tangent vector `X` at identity on the group manifold with the +[`CartanSchoutenZero`](@ref) connection in the direction `d`. See [^Pennec2020] for details. + +[^Pennec2020]: + > X. Pennec and M. Lorenzi, “5 - Beyond Riemannian geometry: The affine connection + > setting for transformation groups,” in Riemannian Geometric Statistics in Medical Image + > Analysis, X. Pennec, S. Sommer, and T. Fletcher, Eds. Academic Press, 2020, pp. 169–229. + > doi: 10.1016/B978-0-12-814725-2.00012-1. +""" +vector_transport_direction( + M::CartanSchoutenZeroGroup, + Y, + ::Identity, + X, + d, + ::ParallelTransport, +) + +function vector_transport_direction!( + M::CartanSchoutenZeroGroup, + Y, + p::Identity, + X, + d, + ::ParallelTransport, +) + dexp_half = group_exp(M.manifold, d / 2) + translate_diff!(M.manifold, Y, dexp_half, p, X, RightAction()) + return translate_diff!(M.manifold, Y, dexp_half, p, Y, LeftAction()) +end + +""" + vector_transport_to(M::CartanSchoutenZeroGroup, ::Identity, X, q, m::ParallelTransport) + +Transport vector `X` at identity of group `M` equipped with the [`CartanSchoutenZero`](@ref) +connection to point `q` using parallel transport. +""" +vector_transport_to(::CartanSchoutenZeroGroup, ::Identity, X, q, ::ParallelTransport) + +function vector_transport_to!( + M::CartanSchoutenZeroGroup, + Y, + p::Identity, + X, + q, + m::ParallelTransport, +) + d = group_log(M.manifold, q) + return vector_transport_direction!(M, Y, p, X, d, m) +end diff --git a/src/groups/group.jl b/src/groups/group.jl index 3f7243e540..d194c0c3b8 100644 --- a/src/groups/group.jl +++ b/src/groups/group.jl @@ -371,6 +371,43 @@ end # Group-specific functions ########################## +@doc raw""" + adjoint_action(G::AbstractGroupManifold, p, X) + +Adjoint action of the element `p` of the Lie group `G` on the element `X` +of the corresponding Lie algebra. + +It is defined as the differential of the group authomorphism ``Ψ_p(q) = pqp⁻¹`` at +the identity of `G`. + +The formula reads +````math +\operatorname{Ad}_p(X) = dΨ_p(e)[X] +```` +where $e$ is the identity element of `G`. + +Note that the adjoint representation of a Lie group isn't generally faithful. +Notably the adjoint representation of SO(2) is trivial. +""" +adjoint_action(G::AbstractGroupManifold, p, X) +@decorator_transparent_function :intransparent function adjoint_action( + G::AbstractGroupManifold, + p, + Xₑ, +) + e = make_identity(G, p) + Xₚ = translate_diff(G, p, e, Xₑ, LeftAction()) + Y = inverse_translate_diff(G, p, p, Xₚ, RightAction()) + return Y +end + +function adjoint_action!(G::AbstractGroupManifold, Y, p, Xₑ) + e = make_identity(G, p) + Xₚ = translate_diff(G, p, e, Xₑ, LeftAction()) + inverse_translate_diff!(G, Y, p, p, Xₚ, RightAction()) + return Y +end + @doc raw""" inv(G::AbstractGroupManifold, p) @@ -441,6 +478,18 @@ end @decorator_transparent_signature compose!(M::AbstractDecoratorManifold, x, p, q) +""" + lie_bracket(G::AbstractGroupManifold, X, Y) + +Lie bracket between elements `X` and `Y` of the Lie algebra corresponding to the Lie group `G`. + +This can be used to compute the adjoint representation of a Lie algebra. +Note that this representation isn't generally faithful. Notably the adjoint +representation of 𝔰𝔬(2) is trivial. +""" +lie_bracket(G::AbstractGroupManifold, X, Y) +@decorator_transparent_signature lie_bracket(M::AbstractDecoratorManifold, X, Y) + _action_order(p, q, ::LeftAction) = (p, q) _action_order(p, q, ::RightAction) = (q, p) @@ -625,7 +674,8 @@ end @doc raw""" group_exp(G::AbstractGroupManifold, X) -Compute the group exponential of the Lie algebra element `X`. +Compute the group exponential of the Lie algebra element `X`. It is equivalent to the +exponential map defined by the [`CartanSchoutenMinus`](@ref) connection. Given an element $X ∈ 𝔤 = T_e \mathcal{G}$, where $e$ is the [`identity`](@ref) element of the group $\mathcal{G}$, and $𝔤$ is its Lie algebra, the group exponential is the map @@ -675,7 +725,8 @@ end @doc raw""" group_log(G::AbstractGroupManifold, q) -Compute the group logarithm of the group element `q`. +Compute the group logarithm of the group element `q`. It is equivalent to the +logarithmic map defined by the [`CartanSchoutenMinus`](@ref) connection. Given an element $q ∈ \mathcal{G}$, compute the right inverse of the group exponential map [`group_exp`](@ref), that is, the element $\log q = X ∈ 𝔤 = T_e \mathcal{G}$, such that @@ -862,6 +913,10 @@ Base.:*(e::Identity{G}, p) where {G<:AdditionGroup} = e Base.:*(p, e::Identity{G}) where {G<:AdditionGroup} = e Base.:*(e::E, ::E) where {G<:AdditionGroup,E<:Identity{G}} = e +adjoint_action(::AdditionGroup, p, X) = X + +adjoint_action!(::AdditionGroup, Y, p, X) = copyto!(Y, X) + Base.zero(e::Identity{G}) where {G<:AdditionGroup} = e Base.identity(::AdditionGroup, p) = zero(p) @@ -899,6 +954,10 @@ group_log(::AdditionGroup, q) = q group_log!(::AdditionGroup, X, q) = copyto!(X, q) +lie_bracket(::AdditionGroup, X, Y) = zero(X) + +lie_bracket!(::AdditionGroup, Z, X, Y) = fill!(Z, 0) + ####################################### # Overloads for MultiplicationOperation ####################################### @@ -974,6 +1033,14 @@ end group_log!(::MultiplicationGroup, X::AbstractMatrix, q::AbstractMatrix) = log_safe!(X, q) +lie_bracket(::MultiplicationGroup, X, Y) = mul!(X * Y, Y, X, -1, true) + +function lie_bracket!(::MultiplicationGroup, Z, X, Y) + mul!(Z, X, Y) + mul!(Z, Y, X, -1, true) + return Z +end + # (a) changes / parent. for f in [ embed, diff --git a/src/groups/semidirect_product_group.jl b/src/groups/semidirect_product_group.jl index 9266cc8f39..2991e767b5 100644 --- a/src/groups/semidirect_product_group.jl +++ b/src/groups/semidirect_product_group.jl @@ -112,6 +112,25 @@ function compose!(G::GT, x, e::E, ::E) where {GT<:SemidirectProductGroup,E<:Iden return identity!(G, x, e) end +@doc raw""" + translate_diff(G::SemidirectProductGroup, p, q, X, conX::LeftAction) + +Perform differential of the left translation on the semidirect product group `G`. + +Since the left translation is defined as (cf. [`SemidirectProductGroup`](@ref)): + +````math +L_{(n', h')} (n, h) = ( L_{n'} θ_{h'}(n), L_{h'} h) +```` + +then its differential can be computed as + +````math +\mathrm{d}L_{(n', h')}(X_n, X_h) = ( \mathrm{d}L_{n'} (\mathrm{d}θ_{h'}(X_n)), \mathrm{d}L_{h'} X_h). +```` +""" +translate_diff(G::SemidirectProductGroup, p, q, X, conX::LeftAction) + function translate_diff!(G::SemidirectProductGroup, Y, p, q, X, conX::LeftAction) M = base_manifold(G) N, H = M.manifolds diff --git a/src/groups/special_euclidean.jl b/src/groups/special_euclidean.jl index 59b52950aa..9c67cfa791 100644 --- a/src/groups/special_euclidean.jl +++ b/src/groups/special_euclidean.jl @@ -84,6 +84,26 @@ Base.@propagate_inbounds function _padvector!( return X end +@doc raw""" + adjoint_action(::SpecialEuclidean{3}, p, fX::TFVector{<:Any,VeeOrthogonalBasis{ℝ}}) + +Adjoint action of the [`SpecialEuclidean`](@ref) group on the vector with coefficients `fX` +tangent at point `p`. + +The formula for the coefficients reads ``t×(R⋅ω) + R⋅r`` for the translation part and +``R⋅ω`` for the rotation part, where `t` is the translation part of `p`, `R` is the rotation +matrix part of `p`, `r` is the translation part of `fX` and `ω` is the rotation part of `fX`, +``×`` is the cross product and ``⋅`` is the matrix product. +""" +function adjoint_action(::SpecialEuclidean{3}, p, fX::TFVector{<:Any,VeeOrthogonalBasis{ℝ}}) + t = p.parts[1] + R = p.parts[2] + r = fX.data[SA[1, 2, 3]] + ω = fX.data[SA[4, 5, 6]] + Rω = R * ω + return TFVector([cross(t, Rω) + R * r; Rω], fX.basis) +end + @doc raw""" affine_matrix(G::SpecialEuclidean, p) -> AbstractMatrix @@ -98,7 +118,14 @@ R & t \\ \end{pmatrix}. ```` +This function embeds $\mathrm{SE}(n)$ in the general linear group $\mathrm{GL}(n+1)$. +It is an isometric embedding and group homomorphism [^RicoMartinez1988]. + See also [`screw_matrix`](@ref) for matrix representations of the Lie algebra. + +[^RicoMartinez1988]: + > Rico Martinez, J. M., “Representations of the Euclidean group and its applications + > to the kinematics of spatial chains,” PhD Thesis, University of Florida, 1988. """ function affine_matrix(G::SpecialEuclidean{n}, p) where {n} pis = submanifold_components(G, p) @@ -128,6 +155,9 @@ the $n + 1 × n + 1$ matrix \end{pmatrix}. ```` +This function embeds $𝔰𝔢(n)$ in the general linear Lie algebra $𝔤𝔩(n+1)$ but it's not +a homomorphic embedding (see [`SpecialEuclideanInGeneralLinear`](@ref) for a homomorphic one). + See also [`affine_matrix`](@ref) for matrix representations of the Lie group. """ function screw_matrix(G::SpecialEuclidean{n}, X) where {n} @@ -386,3 +416,161 @@ function group_log!(G::SpecialEuclidean{3}, X, q) @inbounds _padvector!(G, X) return X end + +""" + lie_bracket(G::SpecialEuclidean, X::ProductRepr, Y::ProductRepr) + lie_bracket(G::SpecialEuclidean, X::AbstractMatrix, Y::AbstractMatrix) + +Calculate the Lie bracket between elements `X` and `Y` of the special Euclidean Lie +algebra. For the matrix representation (which can be obtained using [`screw_matrix`](@ref)) +the formula is ``[X, Y] = XY-YX``, while in the [`ProductRepr`](@ref) representation the +formula reads ``[X, Y] = [(t_1, R_1), (t_2, R_2)] = (R_1 t_2 - R_2 t_1, R_1 R_2 - R_2 R_1)``. +""" +function lie_bracket(G::SpecialEuclidean, X::ProductRepr, Y::ProductRepr) + nX, hX = submanifold_components(G, X) + nY, hY = submanifold_components(G, Y) + return ProductRepr(hX * nY - hY * nX, lie_bracket(G.manifold.manifolds[2], hX, hY)) +end +function lie_bracket(::SpecialEuclidean, X::AbstractMatrix, Y::AbstractMatrix) + return X * Y - Y * X +end + +function lie_bracket!(G::SpecialEuclidean, Z, X, Y) + nX, hX = submanifold_components(G, X) + nY, hY = submanifold_components(G, Y) + nZ, hZ = submanifold_components(G, Z) + lie_bracket!(G.manifold.manifolds[2], hZ, hX, hY) + nZ .= hX * nY .- hY * nX + @inbounds _padvector!(G, Z) + return Z +end + +""" + translate_diff(G::SpecialEuclidean, p, q, X, ::RightAction) + +Differential of the right action of the [`SpecialEuclidean`](@ref) group on itself. +The formula for the rotation part is the differential of the right rotation action, while +the formula for the translation part reads +````math +R_q⋅X_R⋅t_p + X_t +```` +where ``R_q`` is the rotation part of `q`, ``X_R`` is the rotation part of `X`, ``t_p`` +is the translation part of `p` and ``X_t`` is the translation part of `X`. +""" +translate_diff(G::SpecialEuclidean, p, q, X, ::RightAction) + +function translate_diff!(G::SpecialEuclidean, Y, p, q, X, ::RightAction) + np, hp = submanifold_components(G, p) + nq, hq = submanifold_components(G, q) + nX, hX = submanifold_components(G, X) + nY, hY = submanifold_components(G, Y) + hY .= hp' * hX * hp + copyto!(nY, hq * (hX * np) + nX) + @inbounds _padvector!(G, Y) + return Y +end +function translate_diff!(G::SpecialEuclidean, Y, ::Identity, q, X, ::RightAction) + copyto!(G, Y, X) + @inbounds _padvector!(G, Y) + return Y +end +function translate_diff!(G::SpecialEuclidean, Y, p, ::Identity, X, ::RightAction) + np, hp = submanifold_components(G, p) + nX, hX = submanifold_components(G, X) + nY, hY = submanifold_components(G, Y) + hY .= hp' * hX * hp + copyto!(nY, hX * np + nX) + @inbounds _padvector!(G, Y) +end +function translate_diff!(G::SpecialEuclidean, Y, ::Identity, ::Identity, X, ::RightAction) + copyto!(G, Y, X) + @inbounds _padvector!(G, Y) + return Y +end + +@doc raw""" + SpecialEuclideanInGeneralLinear + +An explicit isometric and homomorphic embedding of $\mathrm{SE}(n)$ in $\mathrm{GL}(n+1)$ +and $𝔰𝔢(n)$ in $𝔤𝔩(n+1)$. +Note that this is *not* a transparently isometric embedding. + +# Constructor + + SpecialEuclideanInGeneralLinear(n) +""" +const SpecialEuclideanInGeneralLinear = + EmbeddedManifold{ℝ,<:SpecialEuclidean,<:GeneralLinear} + +function SpecialEuclideanInGeneralLinear(n) + return EmbeddedManifold(SpecialEuclidean(n), GeneralLinear(n + 1)) +end + +""" + embed(M::SpecialEuclideanInGeneralLinear, p) + +Embed the point `p` on [`SpecialEuclidean`](@ref) in the [`GeneralLinear`](@ref) group. +The embedding is calculated using [`affine_matrix`](@ref). +""" +function embed(M::SpecialEuclideanInGeneralLinear, p) + G = M.manifold + return affine_matrix(G, p) +end +""" + embed(M::SpecialEuclideanInGeneralLinear, p, X) + +Embed the tangent vector X at point `p` on [`SpecialEuclidean`](@ref) in the +[`GeneralLinear`](@ref) group. Point `p` can use any representation valid for +`SpecialEuclidean`. The embedding is similar from the one defined by [`screw_matrix`](@ref) +but the translation part is multiplied by inverse of the rotation part. +""" +function embed(M::SpecialEuclideanInGeneralLinear, p, X) + G = M.manifold + np, hp = submanifold_components(G, p) + nX, hX = submanifold_components(G, X) + Y = allocate_result(G, screw_matrix, nX, hX) + nY, hY = submanifold_components(G, Y) + copyto!(hY, hX) + copyto!(nY, hp' * nX) + @inbounds _padvector!(G, Y) + return Y +end + +function embed!(M::SpecialEuclideanInGeneralLinear, q, p) + return copyto!(q, embed(M, p)) +end +function embed!(M::SpecialEuclideanInGeneralLinear, Y, p, X) + return copyto!(Y, embed(M, p, X)) +end + +""" + project(M::SpecialEuclideanInGeneralLinear, p) + +Project point `p` in [`GeneralLinear`](@ref) to the [`SpecialEuclidean`](@ref) group. +This is performed by extracting the rotation and translation part as in [`affine_matrix`](@ref). +""" +function project(M::SpecialEuclideanInGeneralLinear, p) + G = M.manifold + np, hp = submanifold_components(G, p) + return ProductRepr(np, hp) +end +""" + project(M::SpecialEuclideanInGeneralLinear, p, X) + +Project tangent vector `X` at point `p` in [`GeneralLinear`](@ref) to the +[`SpecialEuclidean`](@ref) Lie algebra. +This reverses the transformation performed by [`embed`](@ref embed(M::SpecialEuclideanInGeneralLinear, p, X)) +""" +function project(M::SpecialEuclideanInGeneralLinear, p, X) + G = M.manifold + np, hp = submanifold_components(G, p) + nX, hX = submanifold_components(G, X) + return ProductRepr(hp * nX, hX) +end + +function project!(M::SpecialEuclideanInGeneralLinear, q, p) + return copyto!(q, project(M, p)) +end +function project!(M::SpecialEuclideanInGeneralLinear, Y, p, X) + return copyto!(Y, project(M, p, X)) +end diff --git a/src/manifolds/ConnectionManifold.jl b/src/manifolds/ConnectionManifold.jl new file mode 100644 index 0000000000..2fc85241b2 --- /dev/null +++ b/src/manifolds/ConnectionManifold.jl @@ -0,0 +1,395 @@ +@doc raw""" + AbstractAffineConnection + +Abstract type for affine connections on a manifold. +""" +abstract type AbstractAffineConnection end + +""" + LeviCivitaConnection + +The [Levi-Civita connection](https://en.wikipedia.org/wiki/Levi-Civita_connection) of a Riemannian manifold. +""" +struct LeviCivitaConnection <: AbstractAffineConnection end + +struct MetricDecoratorType <: AbstractDecoratorType end + +""" + AbstractConnectionManifold{𝔽,M<:AbstractManifold{𝔽},G<:AbstractAffineConnection} <: AbstractDecoratorManifold{𝔽} + +Equip an [`AbstractManifold`](@ref) explicitly with an [`AbstractAffineConnection`](@ref) `G`. + +`AbstractConnectionManifold` is defined by values of [`christoffel_symbols_second`](@ref), +which is used for a default implementation of [`exp`](@ref), [`log`](@ref) and +[`vector_transport_to`](@ref). Closed-form formulae for particular connection manifolds may +be explicitly implemented when available. + +An overview of basic properties of affine connection manifolds can be found in [^Pennec2020]. + +[^Pennec2020]: + > X. Pennec and M. Lorenzi, “5 - Beyond Riemannian geometry: The affine connection + > setting for transformation groups,” in Riemannian Geometric Statistics in Medical Image + > Analysis, X. Pennec, S. Sommer, and T. Fletcher, Eds. Academic Press, 2020, pp. 169–229. + > doi: 10.1016/B978-0-12-814725-2.00012-1. +""" +abstract type AbstractConnectionManifold{𝔽} <: + AbstractDecoratorManifold{𝔽,MetricDecoratorType} end + +""" + connection(M::AbstractManifold) + +Get the connection (an object of a subtype of [`AbstractAffineConnection`](@ref)) +of [`AbstractManifold`](@ref) `M`. +""" +connection(::AbstractManifold) + +""" + ConnectionManifold(M, C) + +Decorate the [`AbstractManifold`](@ref) `M` with [`AbstractAffineConnection`](@ref) `C`. +""" +struct ConnectionManifold{𝔽,M<:AbstractManifold{𝔽},C<:AbstractAffineConnection} <: + AbstractConnectionManifold{𝔽} + manifold::M + connection::C +end + +@doc raw""" + christoffel_symbols_second( + M::AbstractManifold, + p, + B::AbstractBasis; + backend::AbstractDiffBackend = diff_backend(), + ) + +Compute the Christoffel symbols of the second kind in local coordinates of basis `B`. +For affine connection manifold the Christoffel symbols need to be explicitly implemented +while, for a [`MetricManifold`](@ref) they are computed as (in Einstein summation convention) + +````math +Γ^{l}_{ij} = g^{kl} Γ_{ijk}, +```` + +where ``Γ_{ijk}`` are the Christoffel symbols of the first kind +(see [`christoffel_symbols_first`](@ref)), and ``g^{kl}`` is the inverse of the local +representation of the metric tensor. The dimensions of the resulting multi-dimensional array +are ordered ``(l,i,j)``. +""" +christoffel_symbols_second(::AbstractManifold, ::Any, ::AbstractBasis) + +@decorator_transparent_signature christoffel_symbols_second( + M::AbstractDecoratorManifold, + p, + B::AbstractBasis; + kwargs..., +) + +@doc raw""" + christoffel_symbols_second_jacobian( + M::AbstractManifold, + p, + B::AbstractBasis; + backend::AbstractDiffBackend = diff_backend(), + ) + +Get partial derivatives of the Christoffel symbols of the second kind +for manifold `M` at `p` with respect to the coordinates of `B`, i.e. + +```math +\frac{∂}{∂ p^l} Γ^{k}_{ij} = Γ^{k}_{ij,l}. +``` + +The dimensions of the resulting multi-dimensional array are ordered ``(i,j,k,l)``. +""" +christoffel_symbols_second_jacobian(::AbstractManifold, ::Any, B::AbstractBasis) +function christoffel_symbols_second_jacobian( + M::AbstractManifold, + p, + B::AbstractBasis; + backend::AbstractDiffBackend=diff_backend(), +) + n = size(p, 1) + ∂Γ = reshape( + _jacobian(q -> christoffel_symbols_second(M, q, B; backend=backend), p, backend), + n, + n, + n, + n, + ) + return ∂Γ +end +@decorator_transparent_signature christoffel_symbols_second_jacobian( + M::AbstractDecoratorManifold, + p, + B::AbstractBasis; + kwargs..., +) + +""" + connection(M::ConnectionManifold) + +Return the connection associated with [`ConnectionManifold`](@ref) `M`. +""" +connection(M::ConnectionManifold) = M.connection + +Base.copyto!(M::AbstractConnectionManifold, q, p) = copyto!(M.manifold, q, p) +Base.copyto!(M::AbstractConnectionManifold, Y, p, X) = copyto!(M.manifold, Y, p, X) + +@doc raw""" + exp(M::AbstractConnectionManifold, p, X) + +Compute the exponential map on the [`AbstractConnectionManifold`](@ref) `M` equipped with +corresponding affine connection. + +If `M` is a [`MetricManifold`](@ref) with a metric that was declared the default metric +using [`is_default_metric`](@ref), this method falls back to `exp(M, p, X)`. + +Otherwise it numerically integrates the underlying ODE, see [`solve_exp_ode`](@ref). +Currently, the numerical integration is only accurate when using a single +coordinate chart that covers the entire manifold. This excludes coordinates +in an embedded space. +""" +exp(::AbstractConnectionManifold, ::Any...) + +@decorator_transparent_fallback function exp!(M::AbstractConnectionManifold, q, p, X) + tspan = (0.0, 1.0) + A = get_default_atlas(M) + i = get_chart_index(M, A, p) + B = induced_basis(M, A, i, TangentSpace) + sol = solve_exp_ode(M, p, X, tspan, B; dense=false, saveat=[1.0]) + n = length(p) + return copyto!(q, sol.u[1][(n + 1):end]) +end + +""" + gaussian_curvature(M::AbstractManifold, p, B::AbstractBasis; backend::AbstractDiffBackend = diff_backend()) + +Compute the Gaussian curvature of the manifold `M` at the point `p` using basis `B`. +This is equal to half of the scalar Ricci curvature, see [`ricci_curvature`](@ref). +""" +gaussian_curvature(::AbstractManifold, ::Any, ::AbstractBasis) +function gaussian_curvature(M::AbstractManifold, p, B::AbstractBasis; kwargs...) + return ricci_curvature(M, p, B; kwargs...) / 2 +end +@decorator_transparent_signature gaussian_curvature( + M::AbstractDecoratorManifold, + p, + B::AbstractBasis; + kwargs..., +) + +function injectivity_radius(M::AbstractConnectionManifold, p) + return injectivity_radius(base_manifold(M), p) +end +function injectivity_radius(M::AbstractConnectionManifold, m::AbstractRetractionMethod) + return injectivity_radius(base_manifold(M), m) +end +function injectivity_radius(M::AbstractConnectionManifold, m::ExponentialRetraction) + return injectivity_radius(base_manifold(M), m) +end +function injectivity_radius(M::AbstractConnectionManifold, p, m::AbstractRetractionMethod) + return injectivity_radius(base_manifold(M), p, m) +end +function injectivity_radius(M::AbstractConnectionManifold, p, m::ExponentialRetraction) + return injectivity_radius(base_manifold(M), p, m) +end + +""" + ricci_tensor(M::AbstractManifold, p, B::AbstractBasis; backend::AbstractDiffBackend = diff_backend()) + +Compute the Ricci tensor, also known as the Ricci curvature tensor, +of the manifold `M` at the point `p` using basis `B`, +see [`https://en.wikipedia.org/wiki/Ricci_curvature#Introduction_and_local_definition`](https://en.wikipedia.org/wiki/Ricci_curvature#Introduction_and_local_definition). +""" +ricci_tensor(::AbstractManifold, ::Any, ::AbstractBasis) +function ricci_tensor(M::AbstractManifold, p, B::AbstractBasis; kwargs...) + R = riemann_tensor(M, p, B; kwargs...) + n = size(R, 1) + Ric = allocate(R, Size(n, n)) + @einsum Ric[i, j] = R[l, i, l, j] + return Ric +end +@decorator_transparent_signature ricci_tensor( + M::AbstractDecoratorManifold, + p, + B::AbstractBasis; + kwargs..., +) + +@doc raw""" + riemann_tensor(M::AbstractManifold, p, B::AbstractBasis; backend::AbstractDiffBackend=diff_backend()) + +Compute the Riemann tensor ``R^l_{ijk}``, also known as the Riemann curvature +tensor, at the point `p` in local coordinates defined by `B`. The dimensions of the +resulting multi-dimensional array are ordered ``(l,i,j,k)``. + +The function uses the coordinate expression involving the second Christoffel symbol, +see [`https://en.wikipedia.org/wiki/Riemann_curvature_tensor#Coordinate_expression`](https://en.wikipedia.org/wiki/Riemann_curvature_tensor#Coordinate_expression) +for details. + +# See also + +[`christoffel_symbols_second`](@ref), [`christoffel_symbols_second_jacobian`](@ref) +""" +riemann_tensor(::AbstractManifold, ::Any, ::AbstractBasis) +function riemann_tensor( + M::AbstractManifold, + p, + B::AbstractBasis; + backend::AbstractDiffBackend=diff_backend(), +) + n = size(p, 1) + Γ = christoffel_symbols_second(M, p, B; backend=backend) + ∂Γ = christoffel_symbols_second_jacobian(M, p, B; backend=backend) ./ n + R = allocate(∂Γ, Size(n, n, n, n)) + @einsum R[l, i, j, k] = + ∂Γ[l, i, k, j] - ∂Γ[l, i, j, k] + Γ[s, i, k] * Γ[l, s, j] - Γ[s, i, j] * Γ[l, s, k] + return R +end +@decorator_transparent_signature riemann_tensor( + M::AbstractDecoratorManifold, + p, + B::AbstractBasis; + kwargs..., +) + +@doc raw""" + solve_exp_ode( + M::AbstractConnectionManifold, + p, + X, + tspan, + B::AbstractBasis; + backend::AbstractDiffBackend = diff_backend(), + solver = AutoVern9(Rodas5()), + kwargs..., + ) + +Approximate the exponential map on the manifold over the provided timespan +assuming the default connection of the given manifold by solving the ordinary differential +equation + +```math +\frac{d^2}{dt^2} p^k + Γ^k_{ij} \frac{d}{dt} p_i \frac{d}{dt} p_j = 0, +``` + +where ``Γ^k_{ij}`` are the Christoffel symbols of the second kind, and +the Einstein summation convention is assumed. The arguments `tspan` and +`solver` follow the `OrdinaryDiffEq` conventions. `kwargs...` specify keyword +arguments that will be passed to `OrdinaryDiffEq.solve`. + +Currently, the numerical integration is only accurate when using a single +coordinate chart that covers the entire manifold. This excludes coordinates +in an embedded space. + +!!! note + This function only works when + [OrdinaryDiffEq.jl](https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl) is loaded with + ```julia + using OrdinaryDiffEq + ``` +""" +function solve_exp_ode(M, p, X, tspan, B::AbstractBasis; kwargs...) + return error( + "solve_exp_ode not implemented on $(typeof(M)) for point $(typeof(p)), vector $(typeof(X)), and timespan $(typeof(tspan)). For a suitable default, enter `using OrdinaryDiffEq` on Julia 1.1 or greater.", + ) +end + +# +# Introduce transparency +# (a) new functions & other parents +for f in [ + christoffel_symbols_second_jacobian, + exp, + gaussian_curvature, + flat, + get_coordinates, + get_vector, + log, + mean, + median, + project, + ricci_tensor, + riemann_tensor, + sharp, + vector_transport_along, + vector_transport_direction, + vector_transport_direction!, #since it has a default using _to! + vector_transport_to, +] + eval( + quote + function decorator_transparent_dispatch( + ::typeof($f), + ::AbstractConnectionManifold, + args..., + ) + return Val(:parent) + end + end, + ) +end + +# (b) changes / intransparencies. +for f in [ + christoffel_symbols_second, # this is basic for connection manifolds but not for metric manifolds + exp!, + flat!, + get_coordinates!, + get_vector!, + get_basis, + inner, + inverse_retract!, + log!, + mean!, + median!, + norm, + project!, + sharp!, + retract!, + vector_transport_along!, + vector_transport_to!, +] + eval( + quote + function decorator_transparent_dispatch( + ::typeof($f), + ::AbstractConnectionManifold, + args..., + ) + return Val(:intransparent) + end + end, + ) +end +# (c) special cases +function decorator_transparent_dispatch( + ::typeof(exp!), + M::AbstractConnectionManifold, + q, + p, + X, + t, +) + return Val(:parent) +end +function decorator_transparent_dispatch( + ::typeof(inverse_retract!), + M::AbstractConnectionManifold, + X, + p, + q, + m::LogarithmicInverseRetraction, +) + return Val(:parent) +end +function decorator_transparent_dispatch( + ::typeof(retract!), + M::AbstractConnectionManifold, + q, + p, + X, + m::ExponentialRetraction, +) + return Val(:parent) +end diff --git a/src/manifolds/MetricManifold.jl b/src/manifolds/MetricManifold.jl index be878bd641..066b00760c 100644 --- a/src/manifolds/MetricManifold.jl +++ b/src/manifolds/MetricManifold.jl @@ -14,8 +14,6 @@ where a zero parameter constructor `T()` is availabe. """ abstract type AbstractMetric end -struct MetricDecoratorType <: AbstractDecoratorType end - # piping syntax for decoration (metric::AbstractMetric)(M::AbstractManifold) = MetricManifold(M, metric) (::Type{T})(M::AbstractManifold) where {T<:AbstractMetric} = MetricManifold(M, T()) @@ -39,7 +37,7 @@ you can of course still implement that directly. Generate the [`AbstractManifold`](@ref) `M` as a manifold with the [`AbstractMetric`](@ref) `G`. """ struct MetricManifold{𝔽,M<:AbstractManifold{𝔽},G<:AbstractMetric} <: - AbstractDecoratorManifold{𝔽,MetricDecoratorType} + AbstractConnectionManifold{𝔽} manifold::M metric::G end @@ -92,26 +90,6 @@ end kwargs..., ) -@doc raw""" - christoffel_symbols_second( - M::AbstractManifold, - p, - B::AbstractBasis; - backend::AbstractDiffBackend = diff_backend(), - ) - -Compute the Christoffel symbols of the second kind in local coordinates of basis `B`. -The Christoffel symbols are (in Einstein summation convention) - -````math -Γ^{l}_{ij} = g^{kl} Γ_{ijk}, -```` - -where ``Γ_{ijk}`` are the Christoffel symbols of the first kind, and -``g^{kl}`` is the inverse of the local representation of the metric tensor. -The dimensions of the resulting multi-dimensional array are ordered ``(l,i,j)``. -""" -christoffel_symbols_second(::AbstractManifold, ::AbstractBasis, ::Any) function christoffel_symbols_second( M::AbstractManifold, p, @@ -124,56 +102,13 @@ function christoffel_symbols_second( @einsum Γ₂[l, i, j] = Ginv[k, l] * Γ₁[i, j, k] return Γ₂ end -@decorator_transparent_signature christoffel_symbols_second( - M::AbstractDecoratorManifold, - p, - B::AbstractBasis; - kwargs..., -) - -@doc raw""" - christoffel_symbols_second_jacobian( - M::AbstractManifold, - p, - B::AbstractBasis; - backend::AbstractDiffBackend = diff_backend(), - ) - -Get partial derivatives of the Christoffel symbols of the second kind -for manifold `M` at `p` with respect to the coordinates of `B`, i.e. - -```math -\frac{∂}{∂ p^l} Γ^{k}_{ij} = Γ^{k}_{ij,l}. -``` -The dimensions of the resulting multi-dimensional array are ordered ``(i,j,k,l)``. """ -christoffel_symbols_second_jacobian(::AbstractManifold, ::Any, B::AbstractBasis) -function christoffel_symbols_second_jacobian( - M::AbstractManifold, - p, - B::AbstractBasis; - backend::AbstractDiffBackend=diff_backend(), -) - n = size(p, 1) - ∂Γ = reshape( - _jacobian(q -> christoffel_symbols_second(M, q, B; backend=backend), p, backend), - n, - n, - n, - n, - ) - return ∂Γ -end -@decorator_transparent_signature christoffel_symbols_second_jacobian( - M::AbstractDecoratorManifold, - p, - B::AbstractBasis; - kwargs..., -) + connection(::MetricManifold) -Base.copyto!(M::MetricManifold, q, p) = copyto!(M.manifold, q, p) -Base.copyto!(M::MetricManifold, Y, p, X) = copyto!(M.manifold, Y, p, X) +Return the [`LeviCivitaConnection`](@ref) for a metric manifold. +""" +connection(::MetricManifold) = LeviCivitaConnection() @doc raw""" det_local_metric(M::AbstractManifold, p, B::AbstractBasis) @@ -218,31 +153,6 @@ end kwargs..., ) -@doc raw""" - exp(N::MetricManifold{M,G}, p, X) - -Copute the exponential map on the [`AbstractManifold`](@ref) `M` equipped with the [`AbstractMetric`](@ref) `G`. - -If the metric was declared the default metric using [`is_default_metric`](@ref), this method -falls back to `exp(M, p, X)`. - -Otherwise it numerically integrates the underlying ODE, see [`solve_exp_ode`](@ref). -Currently, the numerical integration is only accurate when using a single -coordinate chart that covers the entire manifold. This excludes coordinates -in an embedded space. -""" -exp(::MetricManifold, ::Any...) - -@decorator_transparent_fallback function exp!(M::MetricManifold, q, p, X) - tspan = (0.0, 1.0) - A = get_default_atlas(M) - i = get_chart_index(M, A, p) - B = induced_basis(M, A, i, TangentSpace) - sol = solve_exp_ode(M, p, X, tspan, B; dense=false, saveat=[1.0]) - n = length(p) - return copyto!(q, sol.u[1][(n + 1):end]) -end - @doc raw""" flat(N::MetricManifold{M,G}, p, X::FVector{TangentSpaceType}) @@ -268,38 +178,6 @@ flat(::MetricManifold, ::Any...) return ξ end -""" - gaussian_curvature(M::AbstractManifold, p, B::AbstractBasis; backend::AbstractDiffBackend = diff_backend()) - -Compute the Gaussian curvature of the manifold `M` at the point `p` using basis `B`. -""" -gaussian_curvature(::AbstractManifold, ::Any, ::AbstractBasis) -function gaussian_curvature(M::AbstractManifold, p, B::AbstractBasis; kwargs...) - return ricci_curvature(M, p, B; kwargs...) / 2 -end -@decorator_transparent_signature gaussian_curvature( - M::AbstractDecoratorManifold, - p, - B::AbstractBasis; - kwargs..., -) - -function injectivity_radius(M::MetricManifold, p) - return injectivity_radius(base_manifold(M), p) -end -function injectivity_radius(M::MetricManifold, m::AbstractRetractionMethod) - return injectivity_radius(base_manifold(M), m) -end -function injectivity_radius(M::MetricManifold, m::ExponentialRetraction) - return injectivity_radius(base_manifold(M), m) -end -function injectivity_radius(M::MetricManifold, p, m::AbstractRetractionMethod) - return injectivity_radius(base_manifold(M), p, m) -end -function injectivity_radius(M::MetricManifold, p, m::ExponentialRetraction) - return injectivity_radius(base_manifold(M), p, m) -end - @doc raw""" inverse_local_metric(M::AbstractcManifold, p, B::AbstractBasis) @@ -486,10 +364,16 @@ function metric(M::MetricManifold) return M.metric end -""" +@doc raw""" ricci_curvature(M::AbstractManifold, p, B::AbstractBasis; backend::AbstractDiffBackend = diff_backend()) Compute the Ricci scalar curvature of the manifold `M` at the point `p` using basis `B`. +The curvature is computed as the trace of the Ricci curvature tensor with respect to +the metric, that is ``R=g^{ij}R_{ij}`` where ``R`` is the scalar Ricci curvature at `p`, +``g^{ij}`` is the inverse local metric (see [`inverse_local_metric`](@ref)) at `p` and +``R_{ij}`` is the Riccie curvature tensor, see [`ricci_tensor`](@ref). Both the tensor and +inverse local metric are expressed in local coordinates defined by `B`, and the formula +uses the Einstein summation convention. """ ricci_curvature(::AbstractManifold, ::Any, ::AbstractBasis) function ricci_curvature( @@ -509,55 +393,6 @@ end B::AbstractBasis; kwargs..., ) -""" - ricci_tensor(M::AbstractManifold, p, B::AbstractBasis; backend::AbstractDiffBackend = diff_backend()) - -Compute the Ricci tensor, also known as the Ricci curvature tensor, -of the manifold `M` at the point `p` using basis `B`, -see [https://en.wikipedia.org/wiki/Ricci_curvature#Introduction_and_local_definition](https://en.wikipedia.org/wiki/Ricci_curvature#Introduction_and_local_definition). -""" -ricci_tensor(::AbstractManifold, ::Any, ::AbstractBasis) -function ricci_tensor(M::AbstractManifold, p, B::AbstractBasis; kwargs...) - R = riemann_tensor(M, p, B; kwargs...) - n = size(R, 1) - Ric = allocate(R, Size(n, n)) - @einsum Ric[i, j] = R[l, i, l, j] - return Ric -end -@decorator_transparent_signature ricci_tensor( - M::AbstractDecoratorManifold, - p, - B::AbstractBasis; - kwargs..., -) -@doc raw""" - riemann_tensor(M::AbstractManifold, p, B::AbstractBasis; backend::AbstractDiffBackend = diff_backend()) - -Compute the Riemann tensor ``R^l_{ijk}``, also known as the Riemann curvature -tensor, at the point `p`. The dimensions of the resulting multi-dimensional -array are ordered ``(l,i,j,k)``. -""" -riemann_tensor(::AbstractManifold, ::Any, ::AbstractBasis) -function riemann_tensor( - M::AbstractManifold, - p, - B::AbstractBasis; - backend::AbstractDiffBackend=diff_backend(), -) - n = size(p, 1) - Γ = christoffel_symbols_second(M, p, B; backend=backend) - ∂Γ = christoffel_symbols_second_jacobian(M, p, B; backend=backend) ./ n - R = allocate(∂Γ, Size(n, n, n, n)) - @einsum R[l, i, j, k] = - ∂Γ[l, i, k, j] - ∂Γ[l, i, j, k] + Γ[s, i, k] * Γ[l, s, j] - Γ[s, i, j] * Γ[l, s, k] - return R -end -@decorator_transparent_signature riemann_tensor( - M::AbstractDecoratorManifold, - p, - B::AbstractBasis; - kwargs..., -) @doc raw""" sharp(N::MetricManifold{M,G}, p, ξ::FVector{CotangentSpaceType}) @@ -584,84 +419,24 @@ function Base.show(io::IO, M::MetricManifold) return print(io, "MetricManifold($(M.manifold), $(M.metric))") end -@doc raw""" - solve_exp_ode( - M::MetricManifold, - p, - X, - tspan, - B::AbstractBasis; - backend::AbstractDiffBackend = diff_backend(), - solver = AutoVern9(Rodas5()), - kwargs..., - ) - -Approximate the exponential map on the manifold over the provided timespan -assuming the Levi-Civita connection by solving the ordinary differential -equation - -```math -\frac{d^2}{dt^2} p^k + Γ^k_{ij} \frac{d}{dt} p_i \frac{d}{dt} p_j = 0, -``` - -where ``Γ^k_{ij}`` are the Christoffel symbols of the second kind, and -the Einstein summation convention is assumed. The arguments `tspan` and -`solver` follow the `OrdinaryDiffEq` conventions. `kwargs...` specify keyword -arguments that will be passed to `OrdinaryDiffEq.solve`. - -Currently, the numerical integration is only accurate when using a single -coordinate chart that covers the entire manifold. This excludes coordinates -in an embedded space. - -!!! note - This function only works when - [OrdinaryDiffEq.jl](https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl) is loaded with - ```julia - using OrdinaryDiffEq - ``` -""" -function solve_exp_ode(M, p, X, tspan, B::AbstractBasis; kwargs...) - return error( - "solve_exp_ode not implemented on $(typeof(M)) for point $(typeof(p)), vector $(typeof(X)), and timespan $(typeof(tspan)). For a suitable default, enter `using OrdinaryDiffEq` on Julia 1.1 or greater.", - ) -end - # # Introduce transparency # (a) new functions & other parents for f in [ christoffel_symbols_first, - christoffel_symbols_second, - christoffel_symbols_second_jacobian, det_local_metric, einstein_tensor, - exp, - gaussian_curvature, - flat, - get_coordinates, - get_vector, inverse_local_metric, local_metric, local_metric_jacobian, - log, log_local_metric_density, - mean, - median, - project, ricci_curvature, - ricci_tensor, - riemann_tensor, - sharp, - vector_transport_along, - vector_transport_direction, - vector_transport_direction!, #since it has a default using _to! - vector_transport_to, ] eval( quote function decorator_transparent_dispatch( ::typeof($f), - ::MetricManifold, + ::AbstractConnectionManifold, args..., ) return Val(:parent) @@ -670,56 +445,10 @@ for f in [ ) end -# (b) changes / intransparencies. -for f in [ - exp!, - flat!, - get_coordinates!, - get_vector!, - get_basis, - inner, - inverse_retract!, - log!, - mean!, - median!, - norm, - project!, - sharp!, - retract!, - vector_transport_along!, - vector_transport_to!, -] - eval( - quote - function decorator_transparent_dispatch( - ::typeof($f), - ::MetricManifold, - args..., - ) - return Val(:intransparent) - end - end, - ) -end -# (c) special cases -decorator_transparent_dispatch(::typeof(exp!), M::MetricManifold, q, p, X, t) = Val(:parent) function decorator_transparent_dispatch( - ::typeof(inverse_retract!), - M::MetricManifold, - X, - p, - q, - m::LogarithmicInverseRetraction, -) - return Val(:parent) -end -function decorator_transparent_dispatch( - ::typeof(retract!), - M::MetricManifold, - q, - p, - X, - m::ExponentialRetraction, + ::typeof(christoffel_symbols_second), + ::MetricManifold, + args..., ) return Val(:parent) end diff --git a/src/manifolds/Rotations.jl b/src/manifolds/Rotations.jl index 8649cce657..181822adb9 100644 --- a/src/manifolds/Rotations.jl +++ b/src/manifolds/Rotations.jl @@ -720,6 +720,44 @@ Base.show(io::IO, ::Rotations{N}) where {N} = print(io, "Rotations($(N))") Distributions.support(d::NormalRotationDistribution) = MPointSupport(d.manifold) +@doc raw""" + vector_transport_direction(M::Rotations, p, X, d) + +Compute parallel transport of vector `X` tangent at `p` on the [`Rotations`](@ref) +manifold in the direction `d`. The formula, provided in [^Rentmeesters], reads: + +```math +\mathcal P_{q\gets p}X = q^\mathrm{T}p \operatorname{Exp}(d/2) X \operatorname{Exp}(d/2) +``` +where ``q=\exp_p d``. + +The formula simplifies to identity for 2-D rotations. + +[^Rentmeesters]: + > Rentmeesters Q., “A gradient method for geodesic data fitting on some symmetric + > Riemannian manifolds,” in 2011 50th IEEE Conference on Decision and Control and + > European Control Conference, Dec. 2011, pp. 7141–7146. doi: 10.1109/CDC.2011.6161280. +""" +vector_transport_direction(M::Rotations, p, X, d) + +function vector_transport_direction!(M::Rotations, Y, p, X, d, ::ParallelTransport) + expdhalf = exp(d / 2) + q = exp(M, p, d) + return copyto!(Y, transpose(q) * p * expdhalf * X * expdhalf) +end +function vector_transport_direction!(M::Rotations{2}, Y, p, X, d, ::ParallelTransport) + return copyto!(Y, X) +end + +function vector_transport_to!(M::Rotations, Y, p, X, q, ::ParallelTransport) + d = log(M, p, q) + expdhalf = exp(d / 2) + return copyto!(Y, transpose(q) * p * expdhalf * X * expdhalf) +end +function vector_transport_to!(M::Rotations{2}, Y, p, X, q, ::ParallelTransport) + return copyto!(Y, X) +end + @doc raw""" zero_vector(M::Rotations, p) diff --git a/src/tests/tests_group.jl b/src/tests/tests_group.jl index 1d31e19172..bad48187ab 100644 --- a/src/tests/tests_group.jl +++ b/src/tests/tests_group.jl @@ -2,34 +2,39 @@ test_group( G, g_pts::AbstractVector, - v_pts::AbstractVector = [], - ve_pts::AbstractVector = []; + X_pts::AbstractVector = [], + Xe_pts::AbstractVector = []; atol = 1e-10, test_mutating = true, test_group_exp_log = true, test_diff = false, test_invariance = false, + test_lie_bracket=false, + test_adjoint_action=false, diff_convs = [(), (LeftAction(),), (RightAction(),)], ) Tests general properties of the group `G`, given at least three different points elements of it (contained in `g_pts`). -Optionally, specify `test_diff` to test differentials of translation, using `v_pts`, which +Optionally, specify `test_diff` to test differentials of translation, using `X_pts`, which must contain at least one tangent vector at `g_pts[1]`, and the direction conventions specified in `diff_convs`. +`Xe_pts` should contain tangent vectors at identity for testing Lie algebra operations. If the group is equipped with an invariant metric, `test_invariance` indicates that the invariance should be checked for the provided points. """ function test_group( G, g_pts::AbstractVector, - v_pts::AbstractVector=[], - ve_pts::AbstractVector=[]; + X_pts::AbstractVector=[], + Xe_pts::AbstractVector=[]; atol=1e-10, test_mutating=true, test_group_exp_log=true, test_diff=false, test_invariance=false, + test_lie_bracket=false, + test_adjoint_action=false, diff_convs=[(), (LeftAction(),), (RightAction(),)], ) e = make_identity(G, g_pts[1]) @@ -189,7 +194,7 @@ function test_group( end test_diff && Test.@testset "translation differential" begin - X = v_pts[1] + X = X_pts[1] g21 = compose(G, g_pts[2], g_pts[1]) g12 = compose(G, g_pts[1], g_pts[2]) Test.@test isapprox( @@ -283,51 +288,51 @@ function test_group( test_group_exp_log && Test.@testset "group exp/log properties" begin Test.@testset "e = exp(0)" begin - v = group_log(G, identity(G, g_pts[1])) - g = group_exp(G, v) + X = group_log(G, identity(G, g_pts[1])) + g = group_exp(G, X) Test.@test isapprox(G, make_identity(G, g_pts[1]), g; atol=atol) test_mutating && Test.@testset "mutating" begin - v = allocate(ve_pts[1]) - Test.@test group_log!(G, v, identity(G, g_pts[1])) === v + X = allocate(Xe_pts[1]) + Test.@test group_log!(G, X, identity(G, g_pts[1])) === X g = allocate(g_pts[1]) - Test.@test group_exp!(G, g, v) === g + Test.@test group_exp!(G, g, X) === g Test.@test isapprox(G, make_identity(G, g_pts[1]), g; atol=atol) end end - Test.@testset "v = log(exp(v))" begin - for v in ve_pts - g = group_exp(G, v) + Test.@testset "X = log(exp(X))" begin + for X in Xe_pts + g = group_exp(G, X) Test.@test is_point(G, g; atol=atol) - v2 = group_log(G, g) - Test.@test isapprox(G, make_identity(G, g_pts[1]), v2, v; atol=atol) + X2 = group_log(G, g) + Test.@test isapprox(G, make_identity(G, g_pts[1]), X2, X; atol=atol) end test_mutating && Test.@testset "mutating" begin - for v in ve_pts + for X in Xe_pts g = allocate(g_pts[1]) - Test.@test group_exp!(G, g, v) === g + Test.@test group_exp!(G, g, X) === g Test.@test is_point(G, g; atol=atol) - Test.@test isapprox(G, g, group_exp(G, v); atol=atol) - v2 = allocate(v) - Test.@test group_log!(G, v2, g) === v2 - Test.@test isapprox(G, make_identity(G, g_pts[1]), v2, v; atol=atol) + Test.@test isapprox(G, g, group_exp(G, X); atol=atol) + X2 = allocate(X) + Test.@test group_log!(G, X2, g) === X2 + Test.@test isapprox(G, make_identity(G, g_pts[1]), X2, X; atol=atol) end end end Test.@testset "inv(g) = exp(-log(g))" begin g = g_pts[1] - v = group_log(G, g) - ginv = group_exp(G, -v) + X = group_log(G, g) + ginv = group_exp(G, -X) Test.@test isapprox(G, ginv, inv(G, g); atol=atol) end - Test.@testset "exp(sv)∘exp(tv) = exp((s+t)v)" begin - g1 = group_exp(G, 0.2 * ve_pts[1]) - g2 = group_exp(G, 0.3 * ve_pts[1]) - g12 = group_exp(G, 0.5 * ve_pts[1]) + Test.@testset "exp(sX)∘exp(tX) = exp((s+t)X)" begin + g1 = group_exp(G, 0.2 * Xe_pts[1]) + g2 = group_exp(G, 0.3 * Xe_pts[1]) + g12 = group_exp(G, 0.5 * Xe_pts[1]) g1_g2 = compose(G, g1, g2) g2_g1 = compose(G, g2, g1) isapprox(G, g1_g2, g12; atol=atol) @@ -342,17 +347,17 @@ function test_group( y = retract( G, g_pts[1], - v_pts[1], + X_pts[1], Manifolds.GroupExponentialRetraction(conv...), ) Test.@test is_point(G, y; atol=atol) - v2 = inverse_retract( + X2 = inverse_retract( G, g_pts[1], y, Manifolds.GroupLogarithmicInverseRetraction(conv...), ) - Test.@test isapprox(G, g_pts[1], v2, v_pts[1]; atol=atol) + Test.@test isapprox(G, g_pts[1], X2, X_pts[1]; atol=atol) end test_mutating && Test.@testset "mutating" begin @@ -362,19 +367,19 @@ function test_group( G, y, g_pts[1], - v_pts[1], + X_pts[1], Manifolds.GroupExponentialRetraction(conv...), ) === y Test.@test is_point(G, y; atol=atol) - v2 = allocate(v_pts[1]) + X2 = allocate(X_pts[1]) Test.@test inverse_retract!( G, - v2, + X2, g_pts[1], y, Manifolds.GroupLogarithmicInverseRetraction(conv...), - ) === v2 - Test.@test isapprox(G, g_pts[1], v2, v_pts[1]; atol=atol) + ) === X2 + Test.@test isapprox(G, g_pts[1], X2, X_pts[1]; atol=atol) end end end @@ -385,8 +390,8 @@ function test_group( Test.@test has_approx_invariant_metric( G, g_pts[1], - v_pts[1], - v_pts[end], + X_pts[1], + X_pts[end], g_pts, LeftAction(), ) @@ -397,14 +402,68 @@ function test_group( Test.@test has_approx_invariant_metric( G, g_pts[1], - v_pts[1], - v_pts[end], + X_pts[1], + X_pts[end], g_pts, RightAction(), ) end end end + + test_adjoint_action && Test.@testset "Adjoint action" begin + # linearity + X = Xe_pts[1] + Y = Xe_pts[2] + e = identity(G, X) + Test.@test isapprox( + G, + e, + adjoint_action(G, g_pts[2], X + Y), + adjoint_action(G, g_pts[2], X) + adjoint_action(G, g_pts[2], Y), + ) + # inverse property + Test.@test isapprox( + G, + e, + adjoint_action(G, g_pts[2], adjoint_action(G, inv(G, g_pts[2]), X)), + X, + ) + if test_mutating + Z = allocate(X) + adjoint_action!(G, Z, g_pts[2], X) + Test.@test isapprox(G, e, Z, adjoint_action(G, g_pts[2], X)) + end + + # interaction with Lie bracket + if test_lie_bracket + Test.@test isapprox( + G, + e, + adjoint_action(G, g_pts[2], lie_bracket(G, X, Y)), + lie_bracket( + G, + adjoint_action(G, g_pts[2], X), + adjoint_action(G, g_pts[2], Y), + ), + ) + end + end + + test_lie_bracket && Test.@testset "Lie bracket" begin + # anticommutativity + X = X_pts[1] + Y = X_pts[2] + e = identity(G, X) + Test.@test isapprox(G, e, lie_bracket(G, X, Y), -lie_bracket(G, Y, X)) + + if test_mutating + Z = allocate(X) + lie_bracket!(G, Z, X, Y) + Test.@test isapprox(G, e, Z, lie_bracket(G, X, Y)) + end + end + return nothing end @@ -413,7 +472,7 @@ end A::AbstractGroupAction, a_pts::AbstractVector, m_pts::AbstractVector, - v_pts = []; + X_pts = []; atol = 1e-10, atol_ident_compose = 0, test_optimal_alignment = false, @@ -434,7 +493,7 @@ function test_action( A::AbstractGroupAction, a_pts::AbstractVector, m_pts::AbstractVector, - v_pts=[]; + X_pts=[]; atol=1e-10, atol_ident_compose=0, test_optimal_alignment=false, @@ -591,52 +650,52 @@ function test_action( end test_diff && Test.@testset "apply differential" begin - for (m, v) in zip(m_pts, v_pts) + for (m, X) in zip(m_pts, X_pts) for a in a_pts - am, av = apply(A, a, m), apply_diff(A, a, m, v) - ainvm, ainvv = inverse_apply(A, a, m), inverse_apply_diff(A, a, m, v) - Test.@test is_vector(M, am, av, true; atol=atol) + am, aX = apply(A, a, m), apply_diff(A, a, m, X) + ainvm, ainvv = inverse_apply(A, a, m), inverse_apply_diff(A, a, m, X) + Test.@test is_vector(M, am, aX, true; atol=atol) Test.@test is_vector(M, ainvm, ainvv, true; atol=atol) end a12 = compose(A, a_pts[1], a_pts[2]) a2m = apply(A, a_pts[2], m) - a12v = apply_diff(A, a12, m, v) - a2v = apply_diff(A, a_pts[2], m, v) - Test.@test isapprox(M, a2m, apply_diff(A, a_pts[1], a2m, a2v), a12v; atol=atol) + a12X = apply_diff(A, a12, m, X) + a2X = apply_diff(A, a_pts[2], m, X) + Test.@test isapprox(M, a2m, apply_diff(A, a_pts[1], a2m, a2X), a12X; atol=atol) - Test.@test isapprox(M, m, apply_diff(A, e, m, v), v; atol=atol) - Test.@test isapprox(M, m, inverse_apply_diff(A, e, m, v), v; atol=atol) + Test.@test isapprox(M, m, apply_diff(A, e, m, X), X; atol=atol) + Test.@test isapprox(M, m, inverse_apply_diff(A, e, m, X), X; atol=atol) end test_mutating && Test.@testset "mutating" begin - for (m, v) in zip(m_pts, v_pts) + for (m, X) in zip(m_pts, X_pts) for a in a_pts am = apply(A, a, m) - av = allocate(v) - Test.@test apply_diff!(A, av, a, m, v) === av + aX = allocate(X) + Test.@test apply_diff!(A, aX, a, m, X) === aX ainvm = inverse_apply(A, a, m) - ainvv = allocate(v) - Test.@test inverse_apply_diff!(A, ainvv, a, m, v) === ainvv - Test.@test is_vector(M, am, av, true; atol=atol) + ainvv = allocate(X) + Test.@test inverse_apply_diff!(A, ainvv, a, m, X) === ainvv + Test.@test is_vector(M, am, aX, true; atol=atol) Test.@test is_vector(M, ainvm, ainvv, true; atol=atol) end a12 = compose(A, a_pts[1], a_pts[2]) a2m = apply(A, a_pts[2], m) a12m = apply(A, a12, m) - a12v, a2v, a1_a2v = allocate(v), allocate(v), allocate(v) - Test.@test apply_diff!(A, a12v, a12, m, v) === a12v - Test.@test apply_diff!(A, a2v, a_pts[2], m, v) === a2v - Test.@test apply_diff!(A, a1_a2v, a_pts[1], a2m, a2v) === a1_a2v - Test.@test isapprox(M, a12m, a1_a2v, a12v; atol=atol) - - ev = allocate(v) - Test.@test apply_diff!(A, ev, e, m, v) === ev - Test.@test isapprox(G, m, ev, v; atol=atol) - ev = allocate(v) - Test.@test inverse_apply_diff!(A, ev, e, m, v) === ev - Test.@test isapprox(G, m, ev, v; atol=atol) + a12X, a2X, a1_a2X = allocate(X), allocate(X), allocate(X) + Test.@test apply_diff!(A, a12X, a12, m, X) === a12X + Test.@test apply_diff!(A, a2X, a_pts[2], m, X) === a2X + Test.@test apply_diff!(A, a1_a2X, a_pts[1], a2m, a2X) === a1_a2X + Test.@test isapprox(M, a12m, a1_a2X, a12X; atol=atol) + + eX = allocate(X) + Test.@test apply_diff!(A, eX, e, m, X) === eX + Test.@test isapprox(G, m, eX, X; atol=atol) + eX = allocate(X) + Test.@test inverse_apply_diff!(A, eX, e, m, X) === eX + Test.@test isapprox(G, m, eX, X; atol=atol) end end end diff --git a/test/groups/array_manifold.jl b/test/groups/array_manifold.jl index 858c0bc663..2d3f25b36b 100644 --- a/test/groups/array_manifold.jl +++ b/test/groups/array_manifold.jl @@ -71,6 +71,16 @@ include("../utils.jl") group_log!(AG, logp, p2) @test isapprox(G, e, logp.value, group_log(G, p)) + @test lie_bracket(AG, X, X2) isa ValidationTVector + Xlb = allocate(X2) + lie_bracket!(AG, Xlb, X, X2) + @test isapprox(G, e, Xlb.value, lie_bracket(G, X, X2.value)) + + @test adjoint_action(AG, p, X) isa ValidationTVector + Xaa = allocate(X2) + adjoint_action!(AG, Xaa, p, X2) + @test isapprox(G, e, Xaa.value, adjoint_action(G, p, X2.value)) + for conv in (LeftAction(), RightAction()) @test translate(AG, p2, q2, conv) isa ValidationMPoint @test isapprox(G, translate(AG, p2, q2, conv).value, translate(G, p, q, conv)) diff --git a/test/groups/circle_group.jl b/test/groups/circle_group.jl index 62cd6e6254..9eb25ca23d 100644 --- a/test/groups/circle_group.jl +++ b/test/groups/circle_group.jl @@ -35,33 +35,37 @@ using Manifolds: invariant_metric_dispatch, default_metric_dispatch @testset "scalar points" begin pts = [1.0 + 0.0im, 0.0 + 1.0im, (1.0 + 1.0im) / √2] - vpts = [0.0 + 0.5im] + Xpts = [0.0 + 0.5im, 0.0 - 1.5im] @test compose(G, pts[2], pts[1]) ≈ pts[2] * pts[1] - @test translate_diff(G, pts[2], pts[1], vpts[1]) ≈ pts[2] * vpts[1] + @test translate_diff(G, pts[2], pts[1], Xpts[1]) ≈ pts[2] * Xpts[1] test_group( G, pts, - vpts, - vpts; + Xpts, + Xpts; test_diff=true, test_mutating=false, test_invariance=true, + test_lie_bracket=true, + test_adjoint_action=true, ) end @testset "vector points" begin pts = [[1.0 + 0.0im], [0.0 + 1.0im], [(1.0 + 1.0im) / √2]] - vpts = [[0.0 + 0.5im]] + Xpts = [[0.0 + 0.5im], [0.0 - 1.5im]] @test compose(G, pts[2], pts[1]) ≈ pts[2] .* pts[1] - @test translate_diff(G, pts[2], pts[1], vpts[1]) ≈ pts[2] .* vpts[1] + @test translate_diff(G, pts[2], pts[1], Xpts[1]) ≈ pts[2] .* Xpts[1] test_group( G, pts, - vpts, - vpts; + Xpts, + Xpts; test_diff=true, test_mutating=true, test_invariance=true, + test_lie_bracket=true, + test_adjoint_action=true, ) end diff --git a/test/groups/connections.jl b/test/groups/connections.jl new file mode 100644 index 0000000000..0a32c3ecb3 --- /dev/null +++ b/test/groups/connections.jl @@ -0,0 +1,41 @@ +include("../utils.jl") +include("group_utils.jl") + +using Manifolds: connection + +@testset "Cartan-Schouten connections" begin + SO3 = SpecialOrthogonal(3) + SO3minus = ConnectionManifold(SO3, CartanSchoutenMinus()) + SO3plus = ConnectionManifold(SO3, CartanSchoutenPlus()) + SO3zero = ConnectionManifold(SO3, CartanSchoutenZero()) + + e = Matrix{Float64}(I, 3, 3) + p = exp(hat(SO3, Identity(SO3, e), [1.0, 2.0, 3.0])) + q = exp(hat(SO3, Identity(SO3, e), [3.0, 4.0, 1.0])) + X = hat(SO3, Identity(SO3, e), [2.0, 3.0, 4.0]) + SO3e = Identity(SO3, e) + + @testset "connection" begin + @test connection(SO3minus) === CartanSchoutenMinus() + @test connection(SO3plus) === CartanSchoutenPlus() + @test connection(SO3zero) === CartanSchoutenZero() + end + + @testset "log/exp" begin + for CSO3 in [SO3minus, SO3plus, SO3zero] + @test isapprox(SO3, exp(CSO3, p, X), exp(SO3, p, X)) + @test isapprox(SO3, p, log(CSO3, p, q), log(SO3, p, q); atol=1e-6) + end + end + + @testset "Parallel transport" begin + @test isapprox(SO3, q, X, vector_transport_to(SO3minus, SO3e, X, q)) + @test isapprox(SO3, q, q * X / q, vector_transport_to(SO3plus, SO3e, X, q)) + @test isapprox( + SO3, + q, + vector_transport_to(SO3, e, X, q), + vector_transport_to(SO3zero, SO3e, X, q), + ) + end +end diff --git a/test/groups/general_linear.jl b/test/groups/general_linear.jl index 3f0fb79656..1a9e3807d6 100644 --- a/test/groups/general_linear.jl +++ b/test/groups/general_linear.jl @@ -116,7 +116,16 @@ using NLsolve for T in types gpts = convert.(T, pts) vgpts = convert.(T, vpts) - test_group(G, gpts, vgpts, vgpts; test_diff=true, test_invariance=true) + test_group( + G, + gpts, + vgpts, + vgpts; + test_diff=true, + test_invariance=true, + test_lie_bracket=true, + test_adjoint_action=true, + ) test_manifold( G, gpts; diff --git a/test/groups/special_euclidean.jl b/test/groups/special_euclidean.jl index 65fbd247d1..8c4a861462 100644 --- a/test/groups/special_euclidean.jl +++ b/test/groups/special_euclidean.jl @@ -1,6 +1,8 @@ include("../utils.jl") include("group_utils.jl") +using ManifoldsBase: VeeOrthogonalBasis + @testset "Special Euclidean group" begin @testset "SpecialEuclidean($n)" for n in (2, 3, 4) G = SpecialEuclidean(n) @@ -11,18 +13,21 @@ include("group_utils.jl") @test submanifold(G, 1) === TranslationGroup(n) @test submanifold(G, 2) === SpecialOrthogonal(n) Rn = Rotations(n) - x = Matrix(I, n, n) + p = Matrix(I, n, n) if n == 2 t = Vector{Float64}.([1:2, 2:3, 3:4]) ω = [[1.0], [2.0], [1.0]] - tuple_pts = [(ti, exp(Rn, x, hat(Rn, x, ωi))) for (ti, ωi) in zip(t, ω)] - tuple_v = ([-1.0, 2.0], hat(Rn, x, [1.0])) + tuple_pts = [(ti, exp(Rn, p, hat(Rn, p, ωi))) for (ti, ωi) in zip(t, ω)] + tuple_X = [([-1.0, 2.0], hat(Rn, p, [1.0])), ([1.0, -2.0], hat(Rn, p, [0.5]))] elseif n == 3 t = Vector{Float64}.([1:3, 2:4, 4:6]) ω = [[1.0, 2.0, 3.0], [3.0, 2.0, 1.0], [1.0, 3.0, 2.0]] - tuple_pts = [(ti, exp(Rn, x, hat(Rn, x, ωi))) for (ti, ωi) in zip(t, ω)] - tuple_v = ([-1.0, 2.0, 1.0], hat(Rn, x, [1.0, 0.5, -0.5])) + tuple_pts = [(ti, exp(Rn, p, hat(Rn, p, ωi))) for (ti, ωi) in zip(t, ω)] + tuple_X = [ + ([-1.0, 2.0, 1.0], hat(Rn, p, [1.0, 0.5, -0.5])), + ([-2.0, 1.0, 0.5], hat(Rn, p, [-1.0, -0.5, 1.1])), + ] else # n == 4 t = Vector{Float64}.([1:4, 2:5, 3:6]) ω = [ @@ -30,8 +35,11 @@ include("group_utils.jl") [3.0, 2.0, 1.0, 1.0, 2.0, 3.0], [1.0, 3.0, 2.0, 1.0, 2.0, 3.0], ] - tuple_pts = [(ti, exp(Rn, x, hat(Rn, x, ωi))) for (ti, ωi) in zip(t, ω)] - tuple_v = ([-1.0, 2.0, 1.0, 3.0], hat(Rn, x, [1.0, 0.5, -0.5, 0.0, 2.0, 1.0])) + tuple_pts = [(ti, exp(Rn, p, hat(Rn, p, ωi))) for (ti, ωi) in zip(t, ω)] + tuple_X = [ + ([-1.0, 2.0, 1.0, 3.0], hat(Rn, p, [1.0, 0.5, -0.5, 0.0, 2.0, 1.0])), + ([-2.0, 1.5, -1.0, 2.0], hat(Rn, p, [1.0, -0.5, 0.5, 1.0, 0.0, 1.0])), + ] end @testset "product point" begin @@ -39,7 +47,7 @@ include("group_utils.jl") for reshaper in reshapers shape_se = Manifolds.ShapeSpecification(reshaper, M.manifolds...) pts = [Manifolds.prod_point(shape_se, tp...) for tp in tuple_pts] - v_pts = [Manifolds.prod_point(shape_se, tuple_v...)] + X_pts = [Manifolds.prod_point(shape_se, tX...) for tX in tuple_X] g1, g2 = pts[1:2] t1, R1 = g1.parts @@ -50,31 +58,31 @@ include("group_utils.jl") tmp = copy(g1) Manifolds._padpoint!(G, tmp) @test tmp == g1 - tmp = copy(v_pts[1]) + tmp = copy(X_pts[1]) Manifolds._padvector!(G, tmp) - @test tmp == v_pts[1] + @test tmp == X_pts[1] - w = translate_diff(G, pts[1], make_identity(G, pts[1]), v_pts[1]) + w = translate_diff(G, pts[1], make_identity(G, pts[1]), X_pts[1]) w2 = allocate(w) w2.parts[1] .= w.parts[1] w2.parts[2] .= pts[1].parts[2] * w.parts[2] @test screw_matrix(G, w2) ≈ - affine_matrix(G, pts[1]) * screw_matrix(G, v_pts[1]) + affine_matrix(G, pts[1]) * screw_matrix(G, X_pts[1]) test_group( G, pts, - v_pts, - v_pts; + X_pts, + X_pts; test_diff=true, - diff_convs=[(), (LeftAction(),)], + diff_convs=[(), (LeftAction(),), (RightAction(),)], ) end end @testset "product repr" begin pts = [ProductRepr(tp...) for tp in tuple_pts] - v_pts = [ProductRepr(tuple_v...)] + X_pts = [ProductRepr(tX...) for tX in tuple_X] g1, g2 = pts[1:2] t1, R1 = g1.parts @@ -87,50 +95,53 @@ include("group_utils.jl") @test affine_matrix(G, make_identity(G, pts[1])) isa SDiagonal{n,Float64} @test affine_matrix(G, make_identity(G, pts[1])) == SDiagonal{n,Float64}(I) - w = translate_diff(G, pts[1], make_identity(G, pts[1]), v_pts[1]) + w = translate_diff(G, pts[1], make_identity(G, pts[1]), X_pts[1]) w2 = allocate(w) w2.parts[1] .= w.parts[1] w2.parts[2] .= pts[1].parts[2] * w.parts[2] w2mat = screw_matrix(G, w2) - @test w2mat ≈ affine_matrix(G, pts[1]) * screw_matrix(G, v_pts[1]) + @test w2mat ≈ affine_matrix(G, pts[1]) * screw_matrix(G, X_pts[1]) @test screw_matrix(G, w2mat) === w2mat test_group( G, pts, - v_pts, - v_pts; + X_pts, + X_pts; test_diff=true, - diff_convs=[(), (LeftAction(),)], + test_lie_bracket=true, + test_adjoint_action=true, + diff_convs=[(), (LeftAction(),), (RightAction(),)], ) end @testset "affine matrix" begin pts = [affine_matrix(G, ProductRepr(tp...)) for tp in tuple_pts] - v_pts = [screw_matrix(G, ProductRepr(tuple_v...))] + X_pts = [screw_matrix(G, ProductRepr(tX...)) for tX in tuple_X] test_group( G, pts, - v_pts, - v_pts; + X_pts, + X_pts; test_diff=true, - diff_convs=[(), (LeftAction(),)], + test_lie_bracket=true, + diff_convs=[(), (LeftAction(),), (RightAction(),)], ) end @testset "hat/vee" begin shape_se = Manifolds.ShapeSpecification(Manifolds.ArrayReshaper(), M.manifolds...) - x = Manifolds.prod_point(shape_se, tuple_pts[1]...) - V = Manifolds.prod_point(shape_se, tuple_v...) - vexp = [V.parts[1]; vee(Rn, x.parts[2], V.parts[2])] - v = vee(G, x, V) + p = Manifolds.prod_point(shape_se, tuple_pts[1]...) + V = Manifolds.prod_point(shape_se, tuple_X[1]...) + vexp = [V.parts[1]; vee(Rn, p.parts[2], V.parts[2])] + v = vee(G, p, V) @test v ≈ vexp - @test hat(G, x, v) ≈ V + @test hat(G, p, v) ≈ V - v = vee(G, affine_matrix(G, x), screw_matrix(G, V)) + v = vee(G, affine_matrix(G, p), screw_matrix(G, V)) @test v ≈ vexp - @test hat(G, affine_matrix(G, x), v) ≈ screw_matrix(G, V) + @test hat(G, affine_matrix(G, p), v) ≈ screw_matrix(G, V) end end @@ -138,4 +149,74 @@ include("group_utils.jl") @test affine_matrix(G, make_identity(G, ones(12, 12))) isa Diagonal{Float64,Vector{Float64}} @test affine_matrix(G, make_identity(G, ones(12, 12))) == Diagonal(ones(11)) + + @testset "Explicit embedding in GL(n+1)" begin + G = SpecialEuclidean(3) + t = Vector{Float64}.([1:3, 2:4, 4:6]) + ω = [[1.0, 2.0, 3.0], [3.0, 2.0, 1.0], [1.0, 3.0, 2.0]] + p = Matrix(I, 3, 3) + Rn = Rotations(3) + pts = [ProductRepr(ti, exp(Rn, p, hat(Rn, p, ωi))) for (ti, ωi) in zip(t, ω)] + X = ProductRepr([-1.0, 2.0, 1.0], hat(Rn, p, [1.0, 0.5, -0.5])) + q = ProductRepr([0.0, 0.0, 0.0], p) + + GL = GeneralLinear(4) + SEGL = EmbeddedManifold(G, GL) + @test Manifolds.SpecialEuclideanInGeneralLinear(3) === SEGL + pts_gl = [embed(SEGL, pp) for pp in pts] + q_gl = embed(SEGL, q) + X_gl = embed(SEGL, pts_gl[1], X) + + q_gl2 = allocate(q_gl) + embed!(SEGL, q_gl2, q) + @test isapprox(SEGL, q_gl2, q_gl) + + q2 = allocate(q) + project!(SEGL, q2, q_gl) + @test isapprox(G, q, q2) + + @test isapprox(G, pts[1], project(SEGL, pts_gl[1])) + @test isapprox(G, pts[1], X, project(SEGL, pts_gl[1], X_gl)) + + X_gl2 = allocate(X_gl) + embed!(SEGL, X_gl2, pts_gl[1], X) + @test isapprox(SEGL, pts_gl[1], X_gl2, X_gl) + + X2 = allocate(X) + project!(SEGL, X2, pts_gl[1], X_gl) + @test isapprox(G, pts[1], X, X2) + + for conv in [LeftAction(), RightAction()] + tpgl = translate(GL, pts_gl[2], pts_gl[1], conv) + tXgl = translate_diff(GL, pts_gl[2], pts_gl[1], X_gl, conv) + tpse = translate(G, pts[2], pts[1], conv) + tXse = translate_diff(G, pts[2], pts[1], X, conv) + @test isapprox(G, tpse, project(SEGL, tpgl)) + @test isapprox(G, tpse, tXse, project(SEGL, tpgl, tXgl)) + + @test isapprox( + G, + pts_gl[1], + X_gl, + translate_diff(G, Identity(G, pts_gl[1]), pts_gl[1], X_gl, conv), + ) + end + end + + @testset "Adjoint action on 𝔰𝔢(3)" begin + G = SpecialEuclidean(3) + t = Vector{Float64}.([1:3, 2:4, 4:6]) + ω = [[1.0, 2.0, 3.0], [3.0, 2.0, 1.0], [1.0, 3.0, 2.0]] + p = Matrix(I, 3, 3) + Rn = Rotations(3) + pts = [ProductRepr(ti, exp(Rn, p, hat(Rn, p, ωi))) for (ti, ωi) in zip(t, ω)] + X = ProductRepr([-1.0, 2.0, 1.0], hat(Rn, p, [1.0, 0.5, -0.5])) + q = ProductRepr([0.0, 0.0, 0.0], p) + + # adjoint action of SE(3) + fX = TFVector(vee(G, q, X), VeeOrthogonalBasis()) + fXp = adjoint_action(G, pts[1], fX) + fXp2 = adjoint_action(G, pts[1], X) + @test isapprox(G, pts[1], hat(G, pts[1], fXp.data), fXp2) + end end diff --git a/test/groups/special_linear.jl b/test/groups/special_linear.jl index 3e10ed3a2e..d92dcdb531 100644 --- a/test/groups/special_linear.jl +++ b/test/groups/special_linear.jl @@ -100,7 +100,16 @@ using NLsolve for T in types gpts = convert.(T, pts) vgpts = convert.(T, vpts) - test_group(G, gpts, vgpts, vgpts; test_diff=true, test_invariance=true) + test_group( + G, + gpts, + vgpts, + vgpts; + test_diff=true, + test_invariance=true, + test_lie_bracket=true, + test_adjoint_action=true, + ) test_manifold( G, gpts; diff --git a/test/groups/special_orthogonal.jl b/test/groups/special_orthogonal.jl index a115df2e82..744eada676 100644 --- a/test/groups/special_orthogonal.jl +++ b/test/groups/special_orthogonal.jl @@ -46,7 +46,16 @@ include("group_utils.jl") @test translate_diff(G, gpts[2], gpts[1], vgpts[1], LeftAction()) ≈ vgpts[1] @test translate_diff(G, gpts[2], gpts[1], vgpts[1], RightAction()) ≈ transpose(gpts[2]) * vgpts[1] * gpts[2] - test_group(G, gpts, vgpts, vgpts; test_diff=true, test_invariance=true) + test_group( + G, + gpts, + vgpts, + vgpts; + test_diff=true, + test_invariance=true, + test_lie_bracket=true, + test_adjoint_action=true, + ) end @testset "Decorator forwards to group" begin diff --git a/test/groups/translation_group.jl b/test/groups/translation_group.jl index 893a6fd821..65b2b3323f 100644 --- a/test/groups/translation_group.jl +++ b/test/groups/translation_group.jl @@ -17,12 +17,21 @@ include("group_utils.jl") @test base_manifold(G) === Euclidean(2, 3) pts = [reshape(i:(i + 5), (2, 3)) for i in 1:3] - vpts = [reshape(-2:3, (2, 3))] + vpts = [reshape(-2:3, (2, 3)), reshape(-1:4, (2, 3))] for T in types gpts = convert.(T, pts) vgpts = convert.(T, vpts) @test compose(G, gpts[1], gpts[2]) ≈ gpts[1] + gpts[2] - test_group(G, gpts, vgpts, vgpts; test_diff=true, test_invariance=true) + test_group( + G, + gpts, + vgpts, + vgpts; + test_diff=true, + test_invariance=true, + test_lie_bracket=true, + test_adjoint_action=true, + ) end end @@ -34,13 +43,22 @@ include("group_utils.jl") @test base_manifold(G) === Euclidean(2, 3; field=ℂ) pts = [reshape(complex.(i:(i + 5), (i + 1):(i + 6)), (2, 3)) for i in 1:3] - vpts = [reshape(complex.(-2:3, -1:4), (2, 3))] + vpts = [reshape(complex.(-2:3, -1:4), (2, 3)), reshape(complex.(-1:4, 0:5), (2, 3))] for T in types gpts = convert.(T, pts) vgpts = convert.(T, vpts) @test compose(G, gpts[1], gpts[2]) ≈ gpts[1] + gpts[2] @test translate_diff(G, gpts[2], gpts[1], vgpts[1]) ≈ vgpts[1] - test_group(G, gpts, vgpts, vgpts; test_diff=true, test_invariance=true) + test_group( + G, + gpts, + vgpts, + vgpts; + test_diff=true, + test_invariance=true, + test_lie_bracket=true, + test_adjoint_action=true, + ) end end end diff --git a/test/metric.jl b/test/metric.jl index a1e9e02acf..ae7d7bb6d6 100644 --- a/test/metric.jl +++ b/test/metric.jl @@ -1,7 +1,7 @@ using FiniteDifferences, ForwardDiff using LinearAlgebra: I using StatsBase: AbstractWeights, pweights -import Manifolds: mean!, median!, InducedBasis, induced_basis, get_chart_index +import Manifolds: mean!, median!, InducedBasis, induced_basis, get_chart_index, connection include("utils.jl") @@ -199,6 +199,7 @@ end @test TestEuclideanMetric()(E) === M @test TestEuclideanMetric(E) === M + @test connection(M) === Manifolds.LeviCivitaConnection() G = Diagonal(1.0:n) invG = inv(G) diff --git a/test/rotations.jl b/test/rotations.jl index 85c1d31e8e..eebcb13f5c 100644 --- a/test/rotations.jl +++ b/test/rotations.jl @@ -23,17 +23,17 @@ include("utils.jl") @testset "vee/hat" begin M = Manifolds.Rotations(2) - v = [1.23] - x = Matrix{Float64}(I, 2, 2) - V = Manifolds.hat(M, x, v) - @test isa(V, AbstractMatrix) - @test norm(M, x, V) / sqrt(2) ≈ norm(v) - @test Manifolds.vee(M, x, V) == v - - V = project(M, x, randn(2, 2)) - v = Manifolds.vee(M, x, V) - @test isa(v, AbstractVector) - @test Manifolds.hat(M, x, v) == V + Xf = [1.23] + p = Matrix{Float64}(I, 2, 2) + X = Manifolds.hat(M, p, Xf) + @test isa(X, AbstractMatrix) + @test norm(M, p, X) / sqrt(2) ≈ norm(Xf) + @test Manifolds.vee(M, p, X) == Xf + + X = project(M, p, randn(2, 2)) + Xf = Manifolds.vee(M, p, X) + @test isa(Xf, AbstractVector) + @test Manifolds.hat(M, p, Xf) == X end for T in types @@ -46,6 +46,7 @@ include("utils.jl") test_injectivity_radius=false, test_project_tangent=true, test_musical_isomorphisms=true, + test_default_vector_transport=true, retraction_methods=retraction_methods, inverse_retraction_methods=inverse_retraction_methods, point_distributions=[Manifolds.normal_rotation_distribution(M, pts[1], 1.0)], @@ -55,19 +56,19 @@ include("utils.jl") ) @testset "log edge cases" begin - v = Manifolds.hat(M, pts[1], [Float64(π)]) - x = exp(M, pts[1], v) - @test isapprox(x, exp(M, pts[1], log(M, pts[1], x))) + X = Manifolds.hat(M, pts[1], [Float64(π)]) + p = exp(M, pts[1], X) + @test isapprox(p, exp(M, pts[1], log(M, pts[1], p))) end - v = log(M, pts[1], pts[2]) - @test norm(M, pts[1], v) ≈ (angles[2] - angles[1]) * sqrt(2) + X = log(M, pts[1], pts[2]) + @test norm(M, pts[1], X) ≈ (angles[2] - angles[1]) * sqrt(2) # check that exp! does not have a side effect q = allocate(pts[1]) copyto!(M, q, pts[1]) - q2 = exp(M, pts[1], v) - exp!(M, q, q, v) + q2 = exp(M, pts[1], X) + exp!(M, q, q, X) @test norm(q - q2) ≈ 0 v14_polar = inverse_retract(M, pts[1], pts[4], Manifolds.PolarInverseRetraction()) @@ -109,6 +110,7 @@ include("utils.jl") test_injectivity_radius=false, test_musical_isomorphisms=true, test_mutating_rand=true, + test_default_vector_transport=true, retraction_methods=retraction_methods, inverse_retraction_methods=inverse_retraction_methods, point_distributions=[ptd], @@ -120,17 +122,17 @@ include("utils.jl") ) @testset "vee/hat" begin - x = Matrix(1.0I, n, n) - v = randn(manifold_dimension(SOn)) - V = Manifolds.hat(SOn, x, v) - @test isa(V, AbstractMatrix) - @test norm(SOn, x, V) / sqrt(2) ≈ norm(v) - @test Manifolds.vee(SOn, x, V) == v - - V = project(SOn, x, randn(n, n)) - v = Manifolds.vee(SOn, x, V) - @test isa(v, AbstractVector) - @test Manifolds.hat(SOn, x, v) ≈ V + p = Matrix(1.0I, n, n) + Xf = randn(manifold_dimension(SOn)) + X = Manifolds.hat(SOn, p, Xf) + @test isa(X, AbstractMatrix) + @test norm(SOn, p, X) / sqrt(2) ≈ norm(Xf) + @test Manifolds.vee(SOn, p, X) == Xf + + X = project(SOn, p, randn(n, n)) + Xf = Manifolds.vee(SOn, p, X) + @test isa(Xf, AbstractVector) + @test Manifolds.hat(SOn, p, Xf) ≈ X end if n == 4 @@ -147,55 +149,55 @@ include("utils.jl") [0, 0, 10, 0, 0, 1] .* 1e-6, # α ⪆ β ⩰ 0 [0, 0, π / 4, 0, 0, π / 4 - 1e-6], # α ⪆ β > 0 ] - for v in vs - @testset "rotation vector $v" begin - V = Manifolds.hat(SOn, Matrix(1.0I, n, n), v) - x = exp(V) - @test x ≈ exp(SOn, one(x), V) - @test ForwardDiff.derivative(t -> exp(SOn, one(x), t * V), 0) ≈ - V - x2 = exp(log(SOn, one(x), x)) - @test isapprox(x, x2; atol=1e-6) + for Xf in vs + @testset "rotation vector $Xf" begin + X = Manifolds.hat(SOn, Matrix(1.0I, n, n), Xf) + p = exp(X) + @test p ≈ exp(SOn, one(p), X) + @test ForwardDiff.derivative(t -> exp(SOn, one(p), t * X), 0) ≈ + X + p2 = exp(log(SOn, one(p), p)) + @test isapprox(p, p2; atol=1e-6) end end end end - v = Matrix( + X = Matrix( Manifolds.hat(SOn, pts[1], π * normalize(randn(manifold_dimension(SOn)))), ) - x = exp(SOn, pts[1], v) - v2 = log(SOn, pts[1], x) - @test x ≈ exp(SOn, pts[1], v2) + p = exp(SOn, pts[1], X) + v2 = log(SOn, pts[1], p) + @test p ≈ exp(SOn, pts[1], v2) end end @testset "Test AbstractManifold Point and Tangent Vector checks" begin M = Manifolds.Rotations(2) - for x in [1, [2.0 0.0; 0.0 1.0], [1.0 0.5; 0.0 1.0]] - @test_throws DomainError is_point(M, x, true) - @test !is_point(M, x) + for p in [1, [2.0 0.0; 0.0 1.0], [1.0 0.5; 0.0 1.0]] + @test_throws DomainError is_point(M, p, true) + @test !is_point(M, p) end - x = one(zeros(2, 2)) - @test is_point(M, x) - @test is_point(M, x, true) - for v in [1, [0.0 1.0; 0.0 0.0]] - @test_throws DomainError is_vector(M, x, v, true) - @test !is_vector(M, x, v) + p = one(zeros(2, 2)) + @test is_point(M, p) + @test is_point(M, p, true) + for X in [1, [0.0 1.0; 0.0 0.0]] + @test_throws DomainError is_vector(M, p, X, true) + @test !is_vector(M, p, X) end - v = [0.0 1.0; -1.0 0.0] - @test is_vector(M, x, v) - @test is_vector(M, x, v, true) + X = [0.0 1.0; -1.0 0.0] + @test is_vector(M, p, X) + @test is_vector(M, p, X, true) end @testset "Project point" begin M = Manifolds.Rotations(2) - x = Matrix{Float64}(I, 2, 2) - x1 = project(M, x) - @test is_point(M, x1, true) + p = Matrix{Float64}(I, 2, 2) + p1 = project(M, p) + @test is_point(M, p1, true) M = Manifolds.Rotations(3) - x = collect(reshape(1.0:9.0, (3, 3))) - x2 = project(M, x) - @test is_point(M, x2, true) + p = collect(reshape(1.0:9.0, (3, 3))) + p2 = project(M, p) + @test is_point(M, p2, true) rng = MersenneTwister(44) x3 = project(M, randn(rng, 3, 3)) diff --git a/test/runtests.jl b/test/runtests.jl index 80d2c881f1..4c28521780 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -12,7 +12,8 @@ include("utils.jl") if VERSION.prerelease == () && !Sys.iswindows() && VERSION < v"1.6.0" @test length(Test.detect_ambiguities(ManifoldsBase)) <= 66 @test length(Test.detect_ambiguities(Manifolds)) == 0 - @test length(our_base_ambiguities()) <= 24 + # this test takes way too long to perform regularly + # @test length(our_base_ambiguities()) <= 24 else @info "Skipping Ambiguity tests for pre-release versions" end @@ -163,6 +164,7 @@ include("utils.jl") include_test("groups/group_operation_action.jl") include_test("groups/rotation_action.jl") include_test("groups/translation_action.jl") + include_test("groups/connections.jl") include_test("groups/metric.jl") include_test("recipes.jl")