diff --git a/NEWS.md b/NEWS.md index 944bc2b8fa..a242a28386 100644 --- a/NEWS.md +++ b/NEWS.md @@ -5,6 +5,41 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). +## [0.10.0] – 2024-08-24 + +### Changed + +* `Distributions.jl`, `RecursiveArrayTools.jl` and `HybridArrays.jl` were moved to weak dependencies to reduce load time and improve extensibility. +* `translate_diff`, `inv_diff` and thus `apply_diff_group`, are available for all the groups with invariant tangent vector storage. +* `SpecialEuclidean` group now has a different default tangent vector representation, the left-invariant one; to get the old representation pass `vectors=HybridTangentRepresentation()` to the constructor of `SpecialEuclidean`. +* `adjoint_action` takes a direction argument; by default it is `LeftAction`. +* `adjoint_action!` is the necessary method to implement in groups with left-invariant tangent vector representation. +* Fixed a few typos in the doc string of the SPD fixed determinant description. +* Random point on the `MultinomialSymmetricPositiveDefinite` manifold was improved to make it more robust. + +### Added + +* Introduced `exp_inv` and `log_inv` based on `exp_lie` and `log_lie`. They are invariant to the group operation. +* A tutorial about usage of group-related functionality. + +### Removed + +* Deprecated bindings: + * `ExtrinsicEstimation()` (should be replaced with `ExtrinsicEstimation(EfficientEstimator())`), + * `Symplectic` (renamed to `SymplecticMatrices`), + * `SymplecticMatrix` (renamed to `SymplecticElement`). + * `AbstractEstimationMethod` (renamed to `AbstractApproximationMethod`). + * `VectorBundleVectorTransport` (renamed to `FiberBundleProductVectorTransport`). + * `rand` on `SymplecticMatrices` and `SymplecticStiefel` no longer accepts `hamiltonian_norm` as an alias for `σ`. + * `mean!` and `median!` no longer accept `extrinsic_method` (should be replaced with `e = ExtrinsicEstimation(extrinsic_method)`). +* As a result of making `Distributions.jl` and `RecursiveArrayTools.jl` weak dependencies the following symbols are no longer exported from `Manifolds.jl`. Essential functionality is still available but distribution-related features may change in the future without a breaking release. + * `ArrayPartition` (`RecursiveArrayTools.jl` needs to be explicitly imported), + * `ProjectedPointDistribution` (not exported), + * `normal_tvector_distribution` (not exported), + * `projected_distribution` (not exported), + * `uniform_distribution` (not exported). +* Ability to create non-real `SymplecticStiefel` and `SymplecticGrassmann` manifolds; essential functionality was missing so it was removed until a more developed version is developed. + ## [0.9.20] – 2024-06-17 ### Added diff --git a/Project.toml b/Project.toml index 3422c26e17..b563a9d71b 100644 --- a/Project.toml +++ b/Project.toml @@ -1,13 +1,11 @@ name = "Manifolds" uuid = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e" authors = ["Seth Axen ", "Mateusz Baran ", "Ronny Bergmann ", "Antoine Levitt "] -version = "0.9.20" +version = "0.10.0" [deps] -Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Einsum = "b7d42ee7-0b51-5a75-98ca-779d3107e4c0" Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" -HybridArrays = "1baab800-613f-4b0a-84e4-9cd3431bfbb9" Kronecker = "2c470bb0-bcc8-11e8-3dad-c9649493f05e" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" ManifoldDiff = "af67fdf4-a580-4b9f-bbec-742ef357defd" @@ -16,7 +14,6 @@ Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a" MatrixEquations = "99c1a7ee-ab34-5fd5-8076-27c950a045f4" Quaternions = "94ee1d12-ae83-5a48-8b1c-48b8ff168ae0" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" -RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" Requires = "ae029012-a4dd-5104-9daa-d747884805df" SimpleWeightedGraphs = "47aef6b3-ad0c-573a-a1e2-d07658019622" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" @@ -28,17 +25,23 @@ StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" BoundaryValueDiffEq = "764a87c0-6b3e-53db-9096-fe964310641d" Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def" +Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" +HybridArrays = "1baab800-613f-4b0a-84e4-9cd3431bfbb9" NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" +RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [extensions] ManifoldsBoundaryValueDiffEqExt = "BoundaryValueDiffEq" +ManifoldsDistributionsExt = ["Distributions", "RecursiveArrayTools"] +ManifoldsHybridArraysExt = "HybridArrays" ManifoldsNLsolveExt = "NLsolve" -ManifoldsOrdinaryDiffEqDiffEqCallbacksExt = ["DiffEqCallbacks", "OrdinaryDiffEq"] +ManifoldsOrdinaryDiffEqDiffEqCallbacksExt = ["DiffEqCallbacks", "OrdinaryDiffEq", "RecursiveArrayTools"] ManifoldsOrdinaryDiffEqExt = "OrdinaryDiffEq" ManifoldsRecipesBaseExt = ["Colors", "RecipesBase"] +ManifoldsRecursiveArrayToolsExt = "RecursiveArrayTools" ManifoldsTestExt = "Test" [compat] @@ -76,9 +79,11 @@ julia = "1.6" BoundaryValueDiffEq = "764a87c0-6b3e-53db-9096-fe964310641d" Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def" +Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78" FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000" Gtk = "4c0ca9eb-093a-5379-98c5-f87ac0bbbf44" +HybridArrays = "1baab800-613f-4b0a-84e4-9cd3431bfbb9" ImageIO = "82e4d734-157c-48bb-816b-45c225c6df19" ImageMagick = "6218d12a-5da1-5696-b52f-db25d2ecc6d1" NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" @@ -90,8 +95,9 @@ PythonPlot = "274fc56d-3b97-40fa-a1cd-1b4a50311bf9" QuartzImageIO = "dca85d43-d64c-5e67-8c65-017450d5d020" Quaternions = "94ee1d12-ae83-5a48-8b1c-48b8ff168ae0" RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" +RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" VisualRegressionTests = "34922c18-7c2a-561c-bac1-01e79b2c4c92" [targets] -test = ["Test", "BoundaryValueDiffEq", "Colors", "DiffEqCallbacks", "DoubleFloats", "FiniteDifferences", "Gtk", "ImageIO", "ImageMagick", "OrdinaryDiffEq", "NLsolve", "Plots", "PythonPlot", "Quaternions", "QuartzImageIO", "RecipesBase"] +test = ["Test", "BoundaryValueDiffEq", "Colors", "DiffEqCallbacks", "Distributions", "DoubleFloats", "FiniteDifferences", "Gtk", "HybridArrays", "ImageIO", "ImageMagick", "OrdinaryDiffEq", "NLsolve", "Plots", "PythonPlot", "Quaternions", "QuartzImageIO", "RecipesBase", "RecursiveArrayTools"] diff --git a/docs/Project.toml b/docs/Project.toml index c07c63b0ea..3346d0fe7b 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -18,6 +18,7 @@ OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" PythonPlot = "274fc56d-3b97-40fa-a1cd-1b4a50311bf9" RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" +RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" @@ -36,4 +37,5 @@ ManifoldsBase = "0.15.0" OrdinaryDiffEq = "6" Plots = "1" PythonPlot = "1" +RecursiveArrayTools = "3" StaticArrays = "1.0" diff --git a/docs/make.jl b/docs/make.jl index 6e70c54147..8807884d80 100755 --- a/docs/make.jl +++ b/docs/make.jl @@ -24,6 +24,7 @@ if "--quarto" ∈ ARGS tutorials_folder = (@__DIR__) * "/../tutorials" # instantiate the tutorials environment if necessary Pkg.activate(tutorials_folder) + Pkg.develop(PackageSpec(; path=(@__DIR__) * "/../")) Pkg.resolve() Pkg.instantiate() Pkg.activate(@__DIR__) # but return to the docs one before @@ -110,6 +111,7 @@ makedocs(; "perform Hand gesture analysis" => "tutorials/hand-gestures.md", "integrate on manifolds and handle probability densities" => "tutorials/integration.md", "explore curvature without coordinates" => "tutorials/exploring-curvature.md", + "work with groups" => "tutorials/groups.md", ], "Manifolds" => [ "Basic manifolds" => [ diff --git a/docs/src/manifolds/group.md b/docs/src/manifolds/group.md index 5b8f840651..61fe3c0213 100644 --- a/docs/src/manifolds/group.md +++ b/docs/src/manifolds/group.md @@ -43,6 +43,63 @@ Pages = ["groups/GroupManifold.jl"] Order = [:constant, :type, :function] ``` +### Default Representation of Tangent Vectors + +In most groups, the representation of a tangent vector + ``X`` at the point ``p ∈ \mathcal{G}`` is stored + as a vector ``Y ∈ \mathfrak{g}``. +This helps to compute the derivatives of the composition and inverse. + +To explain this, let us assume that *the group consists of matrices* (this is always possible). +The representation of a tangent vector + ``X`` at the point ``p ∈ \mathcal{G}`` is stored + as the vector ``Y ∈ \mathfrak{g}`` given by + +```math +X = pY +``` + +#### Derivative of the Group Composition on the Left + +The derivative of the composition ``pq`` with respect to ``p`` in +the direction ``X``, tangent at ``p`` is given by + +```math +Xq = pYq = pq(q^{-1}Yq) +``` + +We see that with this representation convention, this derivative is just the +adjoint action of ``q^{-1}`` on the vector ``Y``. + +#### Derivative of the Group Composition on the Right + +For the derivative with respect to ``q`` of the composition ``pq`` at a tangent vector ``X`` at ``q`` +stored as ``Y`` with ``X = qY``, we have similarly + +```math +pX = pqY +``` + +With the representation convention above, this derivative is just the identity. + +#### Derivative of the Group Inverse + +Finally, we look at the derivative of the inverse ``p^{-1}`` at a point ``p`` in a tangent direction ``X`` +at ``p`` with ``X = pY``. +The result is a tangent vector at ``p^{-1}`` given by + +```math +-p^{-1}Xp^{-1} = - Yp^{-1} = -p^{-1}(p Y p^{-1}) +``` + +With the representation convention above, this derivative is thus ``-pYp^{-1}``, that is, the opposite of the adjoint action of ``p`` on the vector ``Y``. + +#### Implication for Creating New Groups + +When you create a new group, +defining the adjoint action alone ([`adjoint_action`](@ref)) +automatically defines all the relevant derivatives above. + ### Generic Operations For groups based on an addition operation or a group operation, several default implementations are provided. diff --git a/docs/src/manifolds/power.md b/docs/src/manifolds/power.md index c36c98040c..957f713194 100644 --- a/docs/src/manifolds/power.md +++ b/docs/src/manifolds/power.md @@ -90,7 +90,7 @@ get_component(M, p, 4) The final representation is the [`NestedReplacingPowerRepresentation`](@extref `ManifoldsBase.NestedReplacingPowerRepresentation`). It is similar to the [`NestedPowerRepresentation`](@extref `ManifoldsBase.NestedPowerRepresentation`) but it does not perform in-place operations on the points on the underlying manifold. The example below uses this representation to store points on a power manifold of the [`SpecialEuclidean`](@ref) group in-line in an `Vector` for improved efficiency. When having a mixture of both, i.e. an array structure that is nested (like [´NestedPowerRepresentation](@ref)) in the sense that the elements of the main vector are immutable, then changing the elements can not be done in an in-place way and hence [`NestedReplacingPowerRepresentation`](@extref `ManifoldsBase.NestedReplacingPowerRepresentation`) has to be used. ```@example 4 -using Manifolds, StaticArrays +using Manifolds, StaticArrays, RecursiveArrayTools R2 = Rotations(2) G = SpecialEuclidean(2) diff --git a/docs/src/manifolds/vector_bundle.md b/docs/src/manifolds/vector_bundle.md index 782007896b..a8d3d9e025 100644 --- a/docs/src/manifolds/vector_bundle.md +++ b/docs/src/manifolds/vector_bundle.md @@ -26,7 +26,7 @@ Order = [:constant, :type, :function] The following code defines a point on the tangent bundle of the sphere ``S^2`` and a tangent vector to that point. ```@example tangent-bundle -using Manifolds +using Manifolds, RecursiveArrayTools M = Sphere(2) TB = TangentBundle(M) p = ArrayPartition([1.0, 0.0, 0.0], [0.0, 1.0, 3.0]) diff --git a/docs/src/references.bib b/docs/src/references.bib index c7eb19a168..9683061c59 100644 --- a/docs/src/references.bib +++ b/docs/src/references.bib @@ -195,22 +195,6 @@ @article{LeBrignantPuechmorel:2019 # C # ---------------------------------------------------------------------------------------- -@book{CheegerEbin:2008, - address = {Providence, R.I}, - title = {Comparison {Theorems} in {Riemannian} {Geometry}}, - isbn = {978-0-8218-4417-5}, - publisher = {American Mathematical Society}, - author = {Cheeger, Jeffrey and Ebin, David G.}, - month = may, - year = {2008}, -} -@book{Chikuse:2003, - DOI = {10.1007/978-0-387-21540-2}, - YEAR = {2003}, - AUTHOR = {Yasuko Chikuse}, - PUBLISHER = {Springer New York}, - TITLE = {Statistics on Special Manifolds} -} @inproceedings{ChakrabortyVemuri:2015, AUTHOR = {Rudrasis Chakraborty and Baba C. Vemuri}, BOOKTITLE = {2015 {IEEE} International Conference on Computer Vision (ICCV)}, @@ -229,6 +213,15 @@ @article{ChakrabortyVemuri:2019 TITLE = {Statistics on the Stiefel manifold: Theory and applications}, JOURNAL = {The Annals of Statistics} } +@book{CheegerEbin:2008, + address = {Providence, R.I}, + title = {Comparison {Theorems} in {Riemannian} {Geometry}}, + isbn = {978-0-8218-4417-5}, + publisher = {American Mathematical Society}, + author = {Cheeger, Jeffrey and Ebin, David G.}, + month = may, + year = {2008}, +} @incollection{ChengHoSalehianVemuri:2016, AUTHOR = {Guang Cheng and Jeffrey Ho and Hesamoddin Salehian and Baba C. Vemuri}, BOOKTITLE = {Riemannian Computing in Computer Vision}, @@ -255,6 +248,23 @@ @article{ChevallierLiLuDunson:2022 URL = {https://doi.org/10.48550/arXiv.2009.01983}, YEAR = {2022} } +@book{Chikuse:2003, + DOI = {10.1007/978-0-387-21540-2}, + YEAR = {2003}, + AUTHOR = {Yasuko Chikuse}, + PUBLISHER = {Springer New York}, + TITLE = {Statistics on Special Manifolds} +} +@book{Chirikjian:2012, + edition = {1}, + series = {Applied and {Numerical} {Harmonic} {Analysis}}, + title = {Stochastic {Models}, {Information} {Theory}, and {Lie} {Groups}, {Volume} 2}, + volume = {2}, + url = {https://link.springer.com/book/10.1007/978-0-8176-4944-9}, + publisher = {Birkhäuser Boston, MA}, + author = {Chirikjian, Gregory S.}, + year = {2012}, +} # # D # ---------------------------------------------------------------------------------------- @@ -341,6 +351,21 @@ @article{GallierXu:2002 URL = {http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.35.3205}, YEAR = {2002} } + +@book{GallierQuaintance:2020, + address = {Cham}, + series = {Geometry and {Computing}}, + title = {Differential {Geometry} and {Lie} {Groups}: {A} {Computational} {Perspective}}, + volume = {12}, + isbn = {978-3-030-46039-6 978-3-030-46040-2}, + url = {http://link.springer.com/10.1007/978-3-030-46040-2}, + language = {en}, + publisher = {Springer International Publishing}, + author = {Gallier, Jean and Quaintance, Jocelyn}, + year = {2020}, + doi = {10.1007/978-3-030-46040-2}, +} + @article{GaoSonAbsilStykel:2021, DOI = {10.1137/20m1348522}, YEAR = {2021}, @@ -733,11 +758,16 @@ @article{Sasaki:1958 TITLE = {On the differential geometry of tangent bundles of Riemannian manifolds}, JOURNAL = {Tohoku Math. J.} } -@book{Suhubi:2013, - AUTHOR = {Suhubi, E}, - PUBLISHER = {Academic Press}, - TITLE = {Exterior Analysis: Using Applications of Differential Forms}, - YEAR = {2013} +@misc{SolaDerayAtchuthan:2021, + title = {A micro {Lie} theory for state estimation in robotics}, + url = {http://arxiv.org/abs/1812.01537}, + author = {Solà, Joan and Deray, Jeremie and Atchuthan, Dinesh}, + month = dec, + year = {2021}, + note = {arXiv: 1812.01537}, + eprint={1812.01537}, + archivePrefix={arXiv}, + primaryClass={cs.RO}, } @book{SrivastavaKlassen:2016, AUTHOR = {Anuj Srivastava and Eric P. Klassen}, @@ -747,6 +777,12 @@ @book{SrivastavaKlassen:2016 TITLE = {Functional and Shape Data Analysis}, YEAR = {2016} } +@book{Suhubi:2013, + AUTHOR = {Suhubi, E}, + PUBLISHER = {Academic Press}, + TITLE = {Exterior Analysis: Using Applications of Differential Forms}, + YEAR = {2013} +} # # T # ---------------------------------------------------------------------------------------- diff --git a/ext/ManifoldsDistributionsExt/ManifoldsDistributionsExt.jl b/ext/ManifoldsDistributionsExt/ManifoldsDistributionsExt.jl new file mode 100644 index 0000000000..d068a23e5a --- /dev/null +++ b/ext/ManifoldsDistributionsExt/ManifoldsDistributionsExt.jl @@ -0,0 +1,41 @@ +module ManifoldsDistributionsExt + +if isdefined(Base, :get_extension) + using Manifolds + using Distributions + using Random + using LinearAlgebra + + import Manifolds: + normal_rotation_distribution, + normal_tvector_distribution, + projected_distribution, + uniform_distribution + + using Manifolds: get_iterator, get_parameter, _read, _write + + using RecursiveArrayTools: ArrayPartition +else + # imports need to be relative for Requires.jl-based workflows: + # https://github.com/JuliaArrays/ArrayInterface.jl/pull/387 + using ..Manifolds + using ..Distributions + using ..Random + using ..LinearAlgebra + + import ..Manifolds: + normal_rotation_distribution, + normal_tvector_distribution, + projected_distribution, + uniform_distribution + + using ..Manifolds: get_iterator, get_parameter, _read, _write + + using ..RecursiveArrayTools: ArrayPartition +end + +include("distributions.jl") +include("distributions_for_manifolds.jl") +include("projected_distribution.jl") + +end diff --git a/ext/ManifoldsDistributionsExt/distributions.jl b/ext/ManifoldsDistributionsExt/distributions.jl new file mode 100644 index 0000000000..96afa15699 --- /dev/null +++ b/ext/ManifoldsDistributionsExt/distributions.jl @@ -0,0 +1,126 @@ +""" + FVectorvariate + +Structure that subtypes `VariateForm`, indicating that a single sample +is a vector from a fiber of a vector bundle. +""" +struct FVectorvariate <: VariateForm end + +""" + FVectorSupport(space::AbstractManifold, VectorSpaceFiber) + +Value support for vector bundle fiber-valued distributions (values from a fiber of a vector +bundle at a `point` from the given manifold). +For example used for tangent vector-valued distributions. +""" +struct FVectorSupport{TSpace<:VectorSpaceFiber} <: ValueSupport + space::TSpace +end + +""" + FVectorDistribution{TSpace<:VectorSpaceFiber, T} + +An abstract distribution for vector bundle fiber-valued distributions (values from a fiber +of a vector bundle at point `x` from the given manifold). +For example used for tangent vector-valued distributions. +""" +abstract type FVectorDistribution{TSpace<:VectorSpaceFiber} <: + Distribution{FVectorvariate,FVectorSupport{TSpace}} end + +""" + MPointvariate + +Structure that subtypes `VariateForm`, indicating that a single sample +is a point on a manifold. +""" +struct MPointvariate <: VariateForm end + +""" + MPointSupport(M::AbstractManifold) + +Value support for manifold-valued distributions (values from given +[`AbstractManifold`](@extref `ManifoldsBase.AbstractManifold`) `M`). +""" +struct MPointSupport{TM<:AbstractManifold} <: ValueSupport + manifold::TM +end + +""" + MPointDistribution{TM<:AbstractManifold} + +An abstract distribution for points on manifold of type `TM`. +""" +abstract type MPointDistribution{TM<:AbstractManifold} <: + Distribution{MPointvariate,MPointSupport{TM}} end + +""" + support(d::FVectorDistribution) + +Get the object of type `FVectorSupport` for the distribution `d`. +""" +function Distributions.support(::FVectorDistribution) end + +""" + uniform_distribution(M::Grassmann{<:Any,ℝ}, p) + +Uniform distribution on given (real-valued) [`Grassmann`](@ref) `M`. +Specifically, this is the normalized Haar measure on `M`. +Generated points will be of similar type as `p`. + +The implementation is based on Section 2.5.1 in [Chikuse:2003](@cite); +see also Theorem 2.2.2(iii) in [Chikuse:2003](@cite). +""" +function uniform_distribution(M::Grassmann{<:Any,ℝ}, p) + n, k = get_parameter(M.size) + μ = Distributions.Zeros(n, k) + σ = one(eltype(p)) + Σ1 = Distributions.PDMats.ScalMat(n, σ) + Σ2 = Distributions.PDMats.ScalMat(k, σ) + d = MatrixNormal(μ, Σ1, Σ2) + + return ProjectedPointDistribution(M, d, (M, q, p) -> (q .= svd(p).U), p) +end + +""" + uniform_distribution(M::ProjectiveSpace{<:Any,ℝ}, p) + +Uniform distribution on given [`ProjectiveSpace`](@ref) `M`. Generated points will be of +similar type as `p`. +""" +function uniform_distribution(M::ProjectiveSpace{<:Any,ℝ}, p) + d = Distributions.MvNormal(zero(p), 1.0 * I) + return ProjectedPointDistribution(M, d, project!, p) +end + +""" + uniform_distribution(M::Stiefel{<:Any,ℝ}, p) + +Uniform distribution on given (real-valued) [`Stiefel`](@ref) `M`. +Specifically, this is the normalized Haar and Hausdorff measure on `M`. +Generated points will be of similar type as `p`. + +The implementation is based on Section 2.5.1 in [Chikuse:2003](@cite); +see also Theorem 2.2.1(iii) in [Chikuse:2003](@cite). +""" +function uniform_distribution(M::Stiefel{<:Any,ℝ}, p) + n, k = get_parameter(M.size) + μ = Distributions.Zeros(n, k) + σ = one(eltype(p)) + Σ1 = Distributions.PDMats.ScalMat(n, σ) + Σ2 = Distributions.PDMats.ScalMat(k, σ) + d = MatrixNormal(μ, Σ1, Σ2) + + return ProjectedPointDistribution(M, d, project!, p) +end + +""" + uniform_distribution(M::Sphere{n,ℝ}, p) where {n} + +Uniform distribution on given [`Sphere`](@ref) `M`. Generated points will be of +similar type as `p`. +""" +function uniform_distribution(M::Sphere{<:Any,ℝ}, p) + n = get_parameter(M.size)[1] + d = Distributions.MvNormal(zero(p), 1.0 * I) + return ProjectedPointDistribution(M, d, project!, p) +end diff --git a/ext/ManifoldsDistributionsExt/distributions_for_manifolds.jl b/ext/ManifoldsDistributionsExt/distributions_for_manifolds.jl new file mode 100644 index 0000000000..c9719078cb --- /dev/null +++ b/ext/ManifoldsDistributionsExt/distributions_for_manifolds.jl @@ -0,0 +1,227 @@ + +## product manifold + +""" + ProductPointDistribution(M::ProductManifold, distributions) + +Product distribution on manifold `M`, combined from `distributions`. +""" +struct ProductPointDistribution{ + TM<:ProductManifold, + TD<:(NTuple{N,Distribution} where {N}), +} <: MPointDistribution{TM} + manifold::TM + distributions::TD +end + +function ProductPointDistribution(M::ProductManifold, distributions::MPointDistribution...) + return ProductPointDistribution{typeof(M),typeof(distributions)}(M, distributions) +end +function ProductPointDistribution(distributions::MPointDistribution...) + M = ProductManifold(map(d -> support(d).manifold, distributions)...) + return ProductPointDistribution(M, distributions...) +end + +""" + ProductFVectorDistribution([type::VectorSpaceFiber], [x], distrs...) + +Generates a random vector at point `x` from vector space (a fiber of a tangent +bundle) of type `type` using the product distribution of given distributions. + +Vector space type and `x` can be automatically inferred from distributions `distrs`. +""" +struct ProductFVectorDistribution{ + TSpace<:VectorSpaceFiber{<:Any,<:ProductManifold}, + TD<:(NTuple{N,Distribution} where {N}), +} <: FVectorDistribution{TSpace} + type::TSpace + distributions::TD +end + +function ProductFVectorDistribution(distributions::FVectorDistribution...) + M = ProductManifold(map(d -> support(d).space.manifold, distributions)...) + fiber_type = support(distributions[1]).space.fiber_type + if !all(d -> support(d).space.fiber_type == fiber_type, distributions) + error( + "Not all distributions have support in vector spaces of the same type, which is currently not supported", + ) + end + # Probably worth considering sum spaces in the future? + p = ArrayPartition(map(d -> support(d).space.point, distributions)...) + return ProductFVectorDistribution(Fiber(M, p, fiber_type), distributions) +end + +function Random.rand(rng::AbstractRNG, d::ProductFVectorDistribution) + return ArrayPartition(map(d -> rand(rng, d), d.distributions)...) +end +function Random.rand(rng::AbstractRNG, d::ProductPointDistribution) + return ArrayPartition(map(d -> rand(rng, d), d.distributions)...) +end + +function Distributions._rand!( + rng::AbstractRNG, + d::ProductFVectorDistribution, + X::ArrayPartition, +) + map( + (t1, t2) -> Distributions._rand!(rng, t1, t2), + d.distributions, + submanifold_components(d.type.manifold, X), + ) + return X +end +function Distributions._rand!( + rng::AbstractRNG, + d::ProductPointDistribution, + p::ArrayPartition, +) + map( + (t1, t2) -> Distributions._rand!(rng, t1, t2), + d.distributions, + submanifold_components(d.manifold, p), + ) + return p +end + +Distributions.support(d::ProductPointDistribution) = MPointSupport(d.manifold) +function Distributions.support(tvd::ProductFVectorDistribution) + return FVectorSupport(tvd.type) +end + +function uniform_distribution(M::ProductManifold) + return ProductPointDistribution(M, map(uniform_distribution, M.manifolds)) +end +function uniform_distribution(M::ProductManifold, p) + return ProductPointDistribution( + M, + map(uniform_distribution, M.manifolds, submanifold_components(M, p)), + ) +end + +## power manifold + +""" + PowerPointDistribution(M::AbstractPowerManifold, distribution) + +Power distribution on manifold `M`, based on `distribution`. +""" +struct PowerPointDistribution{TM<:AbstractPowerManifold,TD<:MPointDistribution,TX} <: + MPointDistribution{TM} + manifold::TM + distribution::TD + point::TX +end + +""" + PowerFVectorDistribution([type::VectorSpaceFiber], [x], distr) + +Generates a random vector at a `point` from vector space (a fiber of a tangent +bundle) of type `type` using the power distribution of `distr`. + +Vector space type and `point` can be automatically inferred from distribution `distr`. +""" +struct PowerFVectorDistribution{TSpace<:VectorSpaceFiber,TD<:FVectorDistribution} <: + FVectorDistribution{TSpace} + type::TSpace + distribution::TD +end + +function Random.rand(rng::AbstractRNG, d::PowerFVectorDistribution) + fv = zero_vector(d.type.manifold, d.type.point) + Distributions._rand!(rng, d, fv) + return fv +end +function Random.rand(rng::AbstractRNG, d::PowerPointDistribution) + x = allocate_result(d.manifold, rand, d.point) + Distributions._rand!(rng, d, x) + return x +end + +function Distributions._rand!( + rng::AbstractRNG, + d::PowerFVectorDistribution, + v::AbstractArray, +) + PM = d.type.manifold + rep_size = representation_size(PM.manifold) + for i in get_iterator(d.type.manifold) + copyto!(d.distribution.type.point, _read(PM, rep_size, d.type.point, i)) + Distributions._rand!(rng, d.distribution, _read(PM, rep_size, v, i)) + end + return v +end +function Distributions._rand!(rng::AbstractRNG, d::PowerPointDistribution, x::AbstractArray) + M = d.manifold + rep_size = representation_size(M.manifold) + for i in get_iterator(M) + Distributions._rand!(rng, d.distribution, _write(M, rep_size, x, i)) + end + return x +end + +Distributions.support(tvd::PowerFVectorDistribution) = FVectorSupport(tvd.type) +Distributions.support(d::PowerPointDistribution) = MPointSupport(d.manifold) + +# Rotations + +""" + NormalRotationDistribution(M::Rotations, d::Distribution, x::TResult) + +Distribution that returns a random point on the manifold [`Rotations`](@ref) +`M`. Random point is generated using base distribution `d` and the type +of the result is adjusted to `TResult`. + +See [`normal_rotation_distribution`](@ref) for details. +""" +struct NormalRotationDistribution{TResult,TM<:Rotations,TD<:Distribution} <: + MPointDistribution{TM} + manifold::TM + distr::TD +end + +function NormalRotationDistribution( + M::Rotations, + d::Distribution, + x::TResult, +) where {TResult} + return NormalRotationDistribution{TResult,typeof(M),typeof(d)}(M, d) +end + +function normal_rotation_distribution(M::Rotations, p, σ::Real) + n = get_parameter(M.size)[1] + d = Distributions.MvNormal(zeros(n * n), σ * I) + return NormalRotationDistribution(M, d, p) +end + +function _fix_random_rotation(A::AbstractMatrix) + s = diag(sign.(qr(A).R)) + D = Diagonal(s) + C = qr(A).Q * D + if det(C) < 0 + C[:, [1, 2]] = C[:, [2, 1]] + end + return C +end + +function Random.rand( + rng::AbstractRNG, + d::NormalRotationDistribution{TResult,<:Rotations}, +) where {TResult} + n = get_parameter(d.manifold.size)[1] + return if n == 1 + convert(TResult, ones(1, 1)) + else + A = reshape(rand(rng, d.distr), (n, n)) + convert(TResult, _fix_random_rotation(A)) + end +end + +function Distributions._rand!( + rng::AbstractRNG, + d::NormalRotationDistribution, + x::AbstractArray{<:Real}, +) + return copyto!(x, rand(rng, d)) +end + +Distributions.support(d::NormalRotationDistribution) = MPointSupport(d.manifold) diff --git a/src/projected_distribution.jl b/ext/ManifoldsDistributionsExt/projected_distribution.jl similarity index 100% rename from src/projected_distribution.jl rename to ext/ManifoldsDistributionsExt/projected_distribution.jl diff --git a/ext/ManifoldsHybridArraysExt.jl b/ext/ManifoldsHybridArraysExt.jl new file mode 100644 index 0000000000..ed7b2ada04 --- /dev/null +++ b/ext/ManifoldsHybridArraysExt.jl @@ -0,0 +1,45 @@ +module ManifoldsHybridArraysExt + +if isdefined(Base, :get_extension) + using Manifolds + using ManifoldsBase + + using Manifolds: PowerManifoldMultidimensional + using Manifolds: rep_size_to_colons + + using HybridArrays + + import Manifolds: _read +else + # imports need to be relative for Requires.jl-based workflows: + # https://github.com/JuliaArrays/ArrayInterface.jl/pull/387 + using ..Manifolds + using ..ManifoldsBase + + using ..Manifolds: PowerManifoldMultidimensional + using ..Manifolds: rep_size_to_colons + + using ..HybridArrays + + import ..Manifolds: _read +end + +Base.@propagate_inbounds @inline function _read( + ::PowerManifoldMultidimensional, + rep_size::Tuple, + x::HybridArray, + i::Tuple, +) + return x[rep_size_to_colons(rep_size)..., i...] +end +Base.@propagate_inbounds @inline function _read( + ::PowerManifoldMultidimensional, + rep_size::Tuple{}, + x::HybridArray, + i::NTuple{N,Int}, +) where {N} + # disambiguation + return x[i...] +end + +end diff --git a/ext/ManifoldsOrdinaryDiffEqDiffEqCallbacksExt.jl b/ext/ManifoldsOrdinaryDiffEqDiffEqCallbacksExt.jl index 9229598552..da3971d753 100644 --- a/ext/ManifoldsOrdinaryDiffEqDiffEqCallbacksExt.jl +++ b/ext/ManifoldsOrdinaryDiffEqDiffEqCallbacksExt.jl @@ -13,6 +13,8 @@ if isdefined(Base, :get_extension) using DiffEqCallbacks using OrdinaryDiffEq: OrdinaryDiffEq, SciMLBase, Rodas5, AutoVern9, ODEProblem, solve + + using RecursiveArrayTools: ArrayPartition else # imports need to be relative for Requires.jl-based workflows: # https://github.com/JuliaArrays/ArrayInterface.jl/pull/387 @@ -28,6 +30,8 @@ else using ..DiffEqCallbacks using ..OrdinaryDiffEq: OrdinaryDiffEq, SciMLBase, Rodas5, AutoVern9, ODEProblem, solve + + using ..RecursiveArrayTools: ArrayPartition end """ diff --git a/ext/ManifoldsRecursiveArrayToolsExt/ManifoldsRecursiveArrayToolsExt.jl b/ext/ManifoldsRecursiveArrayToolsExt/ManifoldsRecursiveArrayToolsExt.jl new file mode 100644 index 0000000000..3065773cd7 --- /dev/null +++ b/ext/ManifoldsRecursiveArrayToolsExt/ManifoldsRecursiveArrayToolsExt.jl @@ -0,0 +1,134 @@ +module ManifoldsRecursiveArrayToolsExt + +if isdefined(Base, :get_extension) + using Manifolds + using RecursiveArrayTools: ArrayPartition + using StaticArrays + using LinearAlgebra + + using Base.Iterators: repeated + + using Manifolds: + ActionDirectionAndSide, + ColumnwiseMultiplicationAction, + FiberBundleBasisData, + FiberBundleProductVectorTransport, + LeftColumnwiseSpecialEuclideanAction, + LeftInvariantRepresentation, + PowerManifoldNestedReplacing, + SpecialEuclideanIdentity, + SpecialEuclideanInGeneralLinear, + TangentVectorRepresentation, + TypeParameter + + using Manifolds: + bundle_transport_tangent_direction, bundle_transport_tangent_to, _get_parameter + + import Manifolds: + adjoint_Jacobi_field, + allocate, + allocate_result, + apply, + apply!, + apply_diff, + apply_diff!, + apply_diff_group!, + _common_product_translate_diff, + compose, + _compose, + exp_lie, + get_coordinates, + get_vector, + get_vectors, + hat, + identity_element, + inverse_apply, + inverse_apply_diff, + inverse_translate, + inverse_translate_diff, + isapprox, + jacobi_field, + lie_bracket, + optimal_alignment, + project, + translate, + translate_diff, + vector_representation, + _vector_transport_direction, + _vector_transport_to, + vee +else + # imports need to be relative for Requires.jl-based workflows: + # https://github.com/JuliaArrays/ArrayInterface.jl/pull/387 + using ..Manifolds + using ..RecursiveArrayTools: ArrayPartition + using ..StaticArrays + using ..LinearAlgebra + + using ..Base.Iterators: repeated + + using ..Manifolds: + bundle_transport_tangent_direction, bundle_transport_tangent_to, _get_parameter + + using ..Manifolds: + ActionDirectionAndSide, + ColumnwiseMultiplicationAction, + FiberBundleBasisData, + FiberBundleProductVectorTransport, + LeftColumnwiseSpecialEuclideanAction, + LeftInvariantRepresentation, + PowerManifoldNestedReplacing, + SpecialEuclideanIdentity, + SpecialEuclideanInGeneralLinear, + TangentVectorRepresentation, + TypeParameter + + import ..Manifolds: + adjoint_Jacobi_field, + allocate, + allocate_result, + apply, + apply!, + apply_diff, + apply_diff!, + apply_diff_group!, + _common_product_translate_diff, + compose, + _compose, + exp_lie, + get_coordinates, + get_vector, + get_vectors, + hat, + identity_element, + inverse_apply, + inverse_apply_diff, + inverse_translate, + inverse_translate_diff, + isapprox, + jacobi_field, + lie_bracket, + optimal_alignment, + project, + translate, + translate_diff, + vector_representation, + _vector_transport_direction, + _vector_transport_to, + vee +end + +function allocate( + ::PowerManifoldNestedReplacing, + x::AbstractArray{<:ArrayPartition{T,<:NTuple{N,SArray}}}, +) where {T,N} + return similar(x) +end + +include("bundles_rat.jl") +include("ProductManifoldRAT.jl") +include("ProductGroupRAT.jl") +include("special_euclidean_rat.jl") +include("rotation_translation_action_rat.jl") + +end diff --git a/ext/ManifoldsRecursiveArrayToolsExt/ProductGroupRAT.jl b/ext/ManifoldsRecursiveArrayToolsExt/ProductGroupRAT.jl new file mode 100644 index 0000000000..1bdb013219 --- /dev/null +++ b/ext/ManifoldsRecursiveArrayToolsExt/ProductGroupRAT.jl @@ -0,0 +1,172 @@ + +function allocate_result(G::SemidirectProductGroup, ::typeof(identity_element)) + M = base_manifold(G) + N, H = M.manifolds + np = allocate_result(N, identity_element) + hp = allocate_result(H, identity_element) + return ArrayPartition(np, hp) +end + +function _common_product_translate_diff( + G::ProductGroup, + p, + q, + X::ArrayPartition, + conv::ActionDirectionAndSide, +) + M = G.manifold + return ArrayPartition( + map( + translate_diff, + M.manifolds, + submanifold_components(G, p), + submanifold_components(G, q), + submanifold_components(G, X), + repeated(conv), + )..., + ) +end + +function _compose(M::ProductManifold, p::ArrayPartition, q::ArrayPartition) + return ArrayPartition( + map( + compose, + M.manifolds, + submanifold_components(M, p), + submanifold_components(M, q), + )..., + ) +end + +function Base.exp(M::ProductGroup, p::Identity{ProductOperation}, X::ArrayPartition) + return ArrayPartition( + map( + exp, + M.manifold.manifolds, + submanifold_components(M, p), + submanifold_components(M, X), + )..., + ) +end + +function exp_lie(G::ProductGroup, X) + M = G.manifold + return ArrayPartition(map(exp_lie, M.manifolds, submanifold_components(G, X))...) +end + +Base.@propagate_inbounds function Base.getindex( + p::ArrayPartition, + M::ProductGroup, + i::Union{Integer,Colon,AbstractVector,Val}, +) + return getindex(p, base_manifold(M), i) +end +Base.@propagate_inbounds function Base.getindex( + p::ArrayPartition, + M::SemidirectProductGroup, + i::Union{Integer,Colon,AbstractVector,Val}, +) + return getindex(p, base_manifold(M), i) +end + +function identity_element(G::ProductGroup) + M = G.manifold + return ArrayPartition(map(identity_element, M.manifolds)) +end + +function inverse_translate( + G::ProductGroup, + p::ArrayPartition, + q::ArrayPartition, + conv::ActionDirectionAndSide, +) + M = G.manifold + return ArrayPartition( + map( + inverse_translate, + M.manifolds, + submanifold_components(G, p), + submanifold_components(G, q), + repeated(conv), + )..., + ) +end + +function inverse_translate_diff( + G::ProductGroup, + p::ArrayPartition, + q::ArrayPartition, + X::ArrayPartition, + conv::ActionDirectionAndSide, +) + M = G.manifold + return ArrayPartition( + map( + inverse_translate_diff, + M.manifolds, + submanifold_components(G, p), + submanifold_components(G, q), + submanifold_components(G, X), + repeated(conv), + )..., + ) +end + +function Base.log(M::ProductGroup, p::Identity{ProductOperation}, q::ArrayPartition) + return ArrayPartition( + map( + log, + M.manifold.manifolds, + submanifold_components(M, p), + submanifold_components(M, q), + )..., + ) +end + +Base.@propagate_inbounds function Base.setindex!( + q::ArrayPartition, + p, + M::ProductGroup, + i::Union{Integer,Colon,AbstractVector,Val}, +) + return setindex!(q, p, base_manifold(M), i) +end +Base.@propagate_inbounds function Base.setindex!( + q::ArrayPartition, + p, + M::SemidirectProductGroup, + i::Union{Integer,Colon,AbstractVector,Val}, +) + return setindex!(q, p, base_manifold(M), i) +end + +function translate( + M::ProductGroup, + p::ArrayPartition, + q::ArrayPartition, + conv::ActionDirectionAndSide, +) + return ArrayPartition( + map( + translate, + M.manifold.manifolds, + submanifold_components(M, p), + submanifold_components(M, q), + repeated(conv), + )..., + ) +end + +# these isapprox methods are here just to reduce time-to-first-isapprox +function isapprox(G::ProductGroup, p::ArrayPartition, q::ArrayPartition; kwargs...) + return isapprox(G.manifold, p, q; kwargs...) +end +function isapprox( + G::ProductGroup, + p::ArrayPartition, + X::ArrayPartition, + Y::ArrayPartition; + kwargs..., +) + return isapprox(G.manifold, p, X, Y; kwargs...) +end diff --git a/ext/ManifoldsRecursiveArrayToolsExt/ProductManifoldRAT.jl b/ext/ManifoldsRecursiveArrayToolsExt/ProductManifoldRAT.jl new file mode 100644 index 0000000000..831b589788 --- /dev/null +++ b/ext/ManifoldsRecursiveArrayToolsExt/ProductManifoldRAT.jl @@ -0,0 +1,42 @@ + +function adjoint_Jacobi_field( + M::ProductManifold, + p::ArrayPartition, + q::ArrayPartition, + t, + X::ArrayPartition, + β::Tβ, +) where {Tβ} + return ArrayPartition( + map( + adjoint_Jacobi_field, + M.manifolds, + submanifold_components(M, p), + submanifold_components(M, q), + ntuple(_ -> t, length(M.manifolds)), + submanifold_components(M, X), + ntuple(_ -> β, length(M.manifolds)), + )..., + ) +end + +function jacobi_field( + M::ProductManifold, + p::ArrayPartition, + q::ArrayPartition, + t, + X::ArrayPartition, + β::Tβ, +) where {Tβ} + return ArrayPartition( + map( + jacobi_field, + M.manifolds, + submanifold_components(M, p), + submanifold_components(M, q), + ntuple(_ -> t, length(M.manifolds)), + submanifold_components(M, X), + ntuple(_ -> β, length(M.manifolds)), + )..., + ) +end diff --git a/ext/ManifoldsRecursiveArrayToolsExt/bundles_rat.jl b/ext/ManifoldsRecursiveArrayToolsExt/bundles_rat.jl new file mode 100644 index 0000000000..e575d1ed49 --- /dev/null +++ b/ext/ManifoldsRecursiveArrayToolsExt/bundles_rat.jl @@ -0,0 +1,123 @@ + +@inline function allocate_result(M::FiberBundle, f::TF) where {TF} + p = allocate_result(M.manifold, f) + X = allocate_result(Fiber(M.manifold, p, M.type), f) + return ArrayPartition(p, X) +end + +""" + getindex(p::ArrayPartition, M::FiberBundle, s::Symbol) + p[M::FiberBundle, s] + +Access the element(s) at index `s` of a point `p` on a [`FiberBundle`](@ref) `M` by +using the symbols `:point` and `:vector` or `:fiber` for the base and vector or fiber +component, respectively. +""" +@inline function Base.getindex(p::ArrayPartition, M::FiberBundle, s::Symbol) + (s === :point) && return p.x[1] + (s === :vector || s === :fiber) && return p.x[2] + return throw(DomainError(s, "unknown component $s on $M.")) +end + +function get_vector(M::FiberBundle, p::ArrayPartition, c::AbstractVector, B::AbstractBasis) + n = manifold_dimension(M.manifold) + xp1, xp2 = submanifold_components(M, p) + F = Fiber(M.manifold, xp1, M.type) + return ArrayPartition( + get_vector(M.manifold, xp1, c[1:n], B), + get_vector(F, xp2, c[(n + 1):end], B), + ) +end +function get_vector( + M::FiberBundle, + p::ArrayPartition, + c::AbstractVector, + B::CachedBasis{𝔽,<:AbstractBasis{𝔽},<:FiberBundleBasisData}, +) where {𝔽} + n = manifold_dimension(M.manifold) + xp1, xp2 = submanifold_components(M, p) + F = Fiber(M.manifold, xp1, M.type) + return ArrayPartition( + get_vector(M.manifold, xp1, c[1:n], B.data.base_basis), + get_vector(F, xp2, c[(n + 1):end], B.data.fiber_basis), + ) +end + +function get_vectors( + M::FiberBundle, + p::ArrayPartition, + B::CachedBasis{𝔽,<:AbstractBasis{𝔽},<:FiberBundleBasisData}, +) where {𝔽} + xp1, xp2 = submanifold_components(M, p) + zero_m = zero_vector(M.manifold, xp1) + F = Fiber(M.manifold, xp1, M.type) + zero_f = zero_vector(F, xp1) + vs = typeof(ArrayPartition(zero_m, zero_f))[] + for bv in get_vectors(M.manifold, xp1, B.data.base_basis) + push!(vs, ArrayPartition(bv, zero_f)) + end + for bv in get_vectors(F, xp2, B.data.fiber_basis) + push!(vs, ArrayPartition(zero_m, bv)) + end + return vs +end + +""" + setindex!(p::ArrayPartition, val, M::FiberBundle, s::Symbol) + p[M::VectorBundle, s] = val + +Set the element(s) at index `s` of a point `p` on a [`FiberBundle`](@ref) `M` to `val` by +using the symbols `:point` and `:fiber` or `:vector` for the base and fiber or vector +component, respectively. + +!!! note + + The *content* of element of `p` is replaced, not the element itself. +""" +@inline function Base.setindex!(x::ArrayPartition, val, M::FiberBundle, s::Symbol) + if s === :point + return copyto!(x.x[1], val) + elseif s === :vector || s === :fiber + return copyto!(x.x[2], val) + else + throw(DomainError(s, "unknown component $s on $M.")) + end +end + +function _vector_transport_direction( + M::VectorBundle, + p::ArrayPartition, + X, + d, + m::FiberBundleProductVectorTransport, +) + px, pVx = submanifold_components(M.manifold, p) + VXM, VXF = submanifold_components(M.manifold, X) + dx, dVx = submanifold_components(M.manifold, d) + return ArrayPartition( + vector_transport_direction(M.manifold, px, VXM, dx, m.method_horizontal), + bundle_transport_tangent_direction(M, px, pVx, VXF, dx, m.method_vertical), + ) +end + +function _vector_transport_to( + M::VectorBundle, + p::ArrayPartition, + X, + q, + m::FiberBundleProductVectorTransport, +) + px, pVx = submanifold_components(M.manifold, p) + VXM, VXF = submanifold_components(M.manifold, X) + qx, qVx = submanifold_components(M.manifold, q) + return ArrayPartition( + vector_transport_to(M.manifold, px, VXM, qx, m.method_horizontal), + bundle_transport_tangent_to(M, px, pVx, VXF, qx, m.method_vertical), + ) +end + +@inline function Base.view(x::ArrayPartition, M::FiberBundle, s::Symbol) + (s === :point) && return x.x[1] + (s === :vector || s === :fiber) && return x.x[2] + throw(DomainError(s, "unknown component $s on $M.")) +end diff --git a/ext/ManifoldsRecursiveArrayToolsExt/rotation_translation_action_rat.jl b/ext/ManifoldsRecursiveArrayToolsExt/rotation_translation_action_rat.jl new file mode 100644 index 0000000000..f31b693fc2 --- /dev/null +++ b/ext/ManifoldsRecursiveArrayToolsExt/rotation_translation_action_rat.jl @@ -0,0 +1,174 @@ + +""" + apply(::RotationTranslationActionOnVector{LeftAction}, a::ArrayPartition, p) + +Rotate point `p` by `a.x[2]` and translate it by `a.x[1]`. +""" +function apply(::RotationTranslationActionOnVector{LeftAction}, a::ArrayPartition, p) + return a.x[2] * p + a.x[1] +end +""" + apply(::RotationTranslationActionOnVector{RightAction}, a::ArrayPartition, p) + +Translate point `p` by `-a.x[1]` and rotate it by inverse of `a.x[2]`. +""" +function apply(::RotationTranslationActionOnVector{RightAction}, a::ArrayPartition, p) + return a.x[2] \ (p - a.x[1]) +end +function apply!(::RotationTranslationActionOnVector{LeftAction}, q, a::ArrayPartition, p) + mul!(q, a.x[2], p) + q .+= a.x[1] + return q +end + +""" + inverse_apply(::RotationTranslationActionOnVector{LeftAction}, a::ArrayPartition, p) + +Translate point `p` by `-a.x[1]` and rotate it by inverse of `a.x[2]`. +""" +function inverse_apply( + ::RotationTranslationActionOnVector{LeftAction}, + a::ArrayPartition, + p, +) + return a.x[2] \ (p - a.x[1]) +end +""" + inverse_apply(::RotationTranslationActionOnVector{RightAction}, a::ArrayPartition, p) + +Rotate point `p` by `a.x[2]` and translate it by `a.x[1]`. +""" +function inverse_apply( + ::RotationTranslationActionOnVector{RightAction}, + a::ArrayPartition, + p, +) + return a.x[2] * p + a.x[1] +end + +""" + apply_diff( + ::RotationTranslationActionOnVector{LeftAction}, + a::ArrayPartition, + p, + X, + ) + +Compute differential of `apply` on left [`RotationTranslationActionOnVector`](@ref), +with respect to `p`, i.e. left-multiply vector `X` tangent at `p` by `a.x[2]`. +""" +function apply_diff( + ::RotationTranslationActionOnVector{LeftAction}, + a::ArrayPartition, + p, + X, +) + return a.x[2] * X +end +""" + apply_diff( + ::RotationTranslationActionOnVector{RightAction}, + a::ArrayPartition, + p, + X, + ) + +Compute differential of `apply` on right [`RotationTranslationActionOnVector`](@ref), +with respect to `p`, i.e. left-divide vector `X` tangent at `p` by `a.x[2]`. +""" +function apply_diff( + ::RotationTranslationActionOnVector{RightAction}, + a::ArrayPartition, + p, + X, +) + return a.x[2] \ X +end + +function apply_diff!( + ::RotationTranslationActionOnVector{LeftAction}, + Y, + a::ArrayPartition, + p, + X, +) + mul!(Y, a.x[2], X) + return Y +end +function apply_diff!( + ::RotationTranslationActionOnVector{RightAction}, + Y, + a::ArrayPartition, + p, + X, +) + Y .= a.x[2] \ X + return Y +end + +function apply_diff_group!( + ::RotationTranslationActionOnVector{LeftAction}, + Y, + ::SpecialEuclideanIdentity, + X::ArrayPartition, + p, +) + Y .= X.x[2] * p + return Y +end + +function inverse_apply_diff( + ::RotationTranslationActionOnVector{LeftAction}, + a::ArrayPartition, + p, + X, +) + return a.x[2] \ X +end +function inverse_apply_diff( + ::RotationTranslationActionOnVector{RightAction}, + a::ArrayPartition, + p, + X, +) + return a.x[2] * X +end + +### + +function apply(::LeftColumnwiseSpecialEuclideanAction, a::ArrayPartition, p) + return a.x[2] * p .+ a.x[1] +end + +function apply!(::LeftColumnwiseSpecialEuclideanAction, q, a::ArrayPartition, p) + map((qrow, prow) -> mul!(qrow, a.x[2], prow), eachcol(q), eachcol(p)) + q .+= a.x[1] + return q +end + +function inverse_apply(::LeftColumnwiseSpecialEuclideanAction, a::ArrayPartition, p) + return a.x[2] \ (p .- a.x[1]) +end + +@doc raw""" + optimal_alignment(A::LeftColumnwiseSpecialEuclideanAction, p, q) + +Compute optimal alignment of `p` to `q` under the forward left [`ColumnwiseSpecialEuclideanAction`](@ref). +The algorithm, in sequence, computes optimal translation and optimal rotation. +""" +function optimal_alignment( + A::LeftColumnwiseSpecialEuclideanAction{<:AbstractManifold,<:SpecialEuclidean}, + p, + q, +) + N = _get_parameter(A.SE) + tr_opt = mean(q; dims=1) - mean(p; dims=1) + p_moved = p .+ tr_opt + + Ostar = optimal_alignment( + ColumnwiseMultiplicationAction(A.manifold, SpecialOrthogonal(N)), + p_moved, + q, + ) + return ArrayPartition(tr_opt, Ostar) +end diff --git a/ext/ManifoldsRecursiveArrayToolsExt/special_euclidean_rat.jl b/ext/ManifoldsRecursiveArrayToolsExt/special_euclidean_rat.jl new file mode 100644 index 0000000000..5a84d61a68 --- /dev/null +++ b/ext/ManifoldsRecursiveArrayToolsExt/special_euclidean_rat.jl @@ -0,0 +1,132 @@ + +function lie_bracket(G::SpecialEuclidean, X::ArrayPartition, Y::ArrayPartition) + nX, hX = submanifold_components(G, X) + nY, hY = submanifold_components(G, Y) + return ArrayPartition(hX * nY - hY * nX, lie_bracket(G.manifold.manifolds[2], hX, hY)) +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 ArrayPartition(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) + if vector_representation(M.manifold) isa LeftInvariantRepresentation + return ArrayPartition(nX, hX) + else + return ArrayPartition(hp * nX, hX) + end +end + +### Special methods for better performance of selected operations + +function Base.exp( + M::SpecialEuclidean{T,<:HybridTangentRepresentation}, + p::ArrayPartition, + X::ArrayPartition, +) where {T} + M1, M2 = M.manifold.manifolds + return ArrayPartition( + exp(M1.manifold, p.x[1], X.x[1]), + exp(M2.manifold, p.x[2], X.x[2]), + ) +end +function Base.log( + M::SpecialEuclidean{T,<:HybridTangentRepresentation}, + p::ArrayPartition, + q::ArrayPartition, +) where {T} + M1, M2 = M.manifold.manifolds + return ArrayPartition( + log(M1.manifold, p.x[1], q.x[1]), + log(M2.manifold, p.x[2], q.x[2]), + ) +end +function vee( + M::SpecialEuclidean{T,<:HybridTangentRepresentation}, + p::ArrayPartition, + X::ArrayPartition, +) where {T} + M1, M2 = M.manifold.manifolds + return vcat(vee(M1.manifold, p.x[1], X.x[1]), vee(M2.manifold, p.x[2], X.x[2])) +end +function get_coordinates( + M::SpecialEuclidean{T,<:HybridTangentRepresentation}, + p::ArrayPartition, + X::ArrayPartition, + basis::DefaultOrthogonalBasis, +) where {T} + M1, M2 = M.manifold.manifolds + return vcat( + get_coordinates(M1.manifold, p.x[1], X.x[1], basis), + get_coordinates(M2.manifold, p.x[2], X.x[2], basis), + ) +end +function hat( + M::SpecialEuclidean{TypeParameter{Tuple{2}},<:HybridTangentRepresentation}, + p::ArrayPartition, + c::AbstractVector, +) + M1, M2 = M.manifold.manifolds + return ArrayPartition( + get_vector_orthogonal(M1.manifold, p.x[1], c[SOneTo(2)], ℝ), + get_vector_orthogonal(M2.manifold, p.x[2], c[SA[3]], ℝ), + ) +end +function get_vector( + M::SpecialEuclidean{TypeParameter{Tuple{2}},<:HybridTangentRepresentation}, + p::ArrayPartition, + c::AbstractVector, + basis::DefaultOrthogonalBasis, +) + return ArrayPartition( + get_vector(M.manifold.manifolds[1].manifold, p.x[1], c[SOneTo(2)], basis), + get_vector(M.manifold.manifolds[2].manifold, p.x[2], c[SA[3]], basis), + ) +end + +function hat( + M::SpecialEuclidean{TypeParameter{Tuple{3}},<:HybridTangentRepresentation}, + p::ArrayPartition, + c::AbstractVector, +) + M1, M2 = M.manifold.manifolds + return ArrayPartition( + get_vector_orthogonal(M1.manifold, p.x[1], c[SOneTo(3)], ℝ), + get_vector_orthogonal(M2.manifold, p.x[2], c[SA[4, 5, 6]], ℝ), + ) +end +function get_vector( + M::SpecialEuclidean{TypeParameter{Tuple{3}},<:HybridTangentRepresentation}, + p::ArrayPartition, + c::AbstractVector, + basis::DefaultOrthogonalBasis, +) + return ArrayPartition( + get_vector(M.manifold.manifolds[1].manifold, p.x[1], c[SOneTo(3)], basis), + get_vector(M.manifold.manifolds[2].manifold, p.x[2], c[SA[4, 5, 6]], basis), + ) +end +function compose( + ::SpecialEuclidean{T,<:HybridTangentRepresentation}, + p::ArrayPartition, + q::ArrayPartition, +) where {T} + return ArrayPartition(p.x[2] * q.x[1] + p.x[1], p.x[2] * q.x[2]) +end diff --git a/ext/ManifoldsTestExt/tests_general.jl b/ext/ManifoldsTestExt/tests_general.jl index f0cad921fd..cc056bff49 100644 --- a/ext/ManifoldsTestExt/tests_general.jl +++ b/ext/ManifoldsTestExt/tests_general.jl @@ -755,7 +755,6 @@ function test_manifold( for p in pts prand = allocate(p) for pd in point_distributions - Test.@test Manifolds.support(pd) isa Manifolds.MPointSupport{typeof(M)} for _ in 1:10 Test.@test is_point(M, rand(pd)) if test_mutating_rand @@ -829,13 +828,16 @@ function test_manifold( Test.@testset "tangent vector distributions" begin # COV_EXCL_LINE for tvd in tvector_distributions - supp = Manifolds.support(tvd) - Test.@test supp isa - Manifolds.FVectorSupport{<:TangentSpace{number_system(M),typeof(M)}} - for _ in 1:10 + p = tvd.type.point + for _ in 1:5 randtv = rand(tvd) atol = rand_tvector_atol_multiplier * find_eps(randtv) - Test.@test is_vector(M, supp.space.point, randtv, true; atol=atol) + Test.@test is_vector(M, p, randtv, true; atol=atol) + if test_mutating_rand + X = allocate(randtv) + rand!(tvd, X) + Test.@test is_vector(M, p, X, true; atol=atol) + end end end end diff --git a/ext/ManifoldsTestExt/tests_group.jl b/ext/ManifoldsTestExt/tests_group.jl index 1b6023ff44..254f08f1c4 100644 --- a/ext/ManifoldsTestExt/tests_group.jl +++ b/ext/ManifoldsTestExt/tests_group.jl @@ -114,7 +114,13 @@ function test_group( test_adjoint_action::Bool=false, test_inv_diff::Bool=false, test_adjoint_inv_diff::Bool=false, - diff_convs=[(), (LeftForwardAction(),), (RightBackwardAction(),)], + diff_convs=[ + (), + (LeftForwardAction(),), + (RightBackwardAction(),), + (LeftBackwardAction(),), + (RightForwardAction(),), + ], test_log_from_identity::Bool=false, test_exp_from_identity::Bool=false, test_vee_hat_from_identity::Bool=false, @@ -379,6 +385,12 @@ function test_group( Test.@testset "test_inv_diff" for side in [LeftSide(), RightSide()] test_inv_diff_fn(G, g_pts[1], X_pts[1], side) end # COV_EXCL_LINE + if test_mutating + Y = inv_diff(G, e, Xe_pts[1]) + Z = allocate(Y) + inv_diff!(G, Z, e, Xe_pts[1]) + Test.@test isapprox(TangentSpace(G, e), Y, Z) + end end test_adjoint_inv_diff && Test.@testset "Differential of inverse" begin # COV_EXCL_LINE Test.@test isapprox(adjoint_inv_diff(G, e, Xe_pts[1]), -Xe_pts[1]; atol=atol) @@ -534,10 +546,27 @@ function test_group( adjoint_action(G, g_pts[2], adjoint_action(G, inv(G, g_pts[2]), X)), X, ) + # right adjoint action + Test.@test isapprox( + G, + e, + adjoint_action(G, g_pts[2], adjoint_action(G, g_pts[2], X, RightAction())), + X, + ) + # adjoint action at identity + for conv in [LeftAction(), RightAction()] + isapprox( + TangentSpace(G, identity_element(G)), + adjoint_action(G, e, Xe_pts[1], conv), + Xe_pts[1], + ) + end 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)) + adjoint_action!(G, Z, g_pts[2], X, RightAction()) + Test.@test isapprox(G, e, Z, adjoint_action(G, g_pts[2], X, RightAction())) end # interaction with Lie bracket diff --git a/src/Manifolds.jl b/src/Manifolds.jl index 3ae857b7a2..e693a5b882 100644 --- a/src/Manifolds.jl +++ b/src/Manifolds.jl @@ -190,9 +190,7 @@ import Statistics: mean, mean!, median, median!, cov, var import StatsBase: mean_and_var using Base.Iterators: repeated -using Distributions using Einsum: @einsum -using HybridArrays using Kronecker using Graphs using LinearAlgebra @@ -356,7 +354,6 @@ using Markdown: @doc_str using MatrixEquations: lyapc, sylvc using Quaternions: Quaternions using Random -using RecursiveArrayTools: ArrayPartition using Requires using SimpleWeightedGraphs: AbstractSimpleWeightedGraph, get_weight using SpecialFunctions @@ -367,10 +364,28 @@ using StatsBase: AbstractWeights include("utils.jl") -include("product_representations.jl") - include("manifold_fallbacks.jl") +""" + projected_distribution(M::AbstractManifold, d, [p=rand(d)]) + +Wrap the standard distribution `d` into a manifold-valued distribution. Generated +points will be of similar type to `p`. By default, the type is not changed. +""" +function projected_distribution end + +""" + normal_tvector_distribution(M::AbstractManifold, p, σ) + +Normal distribution in ambient space with standard deviation `σ` +projected to tangent space at `p`. +""" +function normal_tvector_distribution end + +function uniform_distribution(M::AbstractManifold) + return uniform_distribution(M, allocate_result(M, uniform_distribution)) +end + # Main Meta Manifolds include("manifolds/ConnectionManifold.jl") include("manifolds/MetricManifold.jl") @@ -382,8 +397,6 @@ include("manifolds/VectorBundle.jl") include("groups/group.jl") # Features I: Which are extended on Meta Manifolds -include("distributions.jl") -include("projected_distribution.jl") include("statistics.jl") # Meta Manifolds II: Products @@ -589,12 +602,6 @@ function __init__() end @static if !isdefined(Base, :get_extension) - @require OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" begin - @require DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def" begin - include("../ext/ManifoldsOrdinaryDiffEqDiffEqCallbacksExt.jl") - end - end - @require OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" begin include("../ext/ManifoldsOrdinaryDiffEqExt.jl") end @@ -616,6 +623,26 @@ function __init__() @require Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" begin include("../ext/ManifoldsTestExt/ManifoldsTestExt.jl") end + + @require RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" begin + include( + "../ext/ManifoldsRecursiveArrayToolsExt/ManifoldsRecursiveArrayToolsExt.jl", + ) + + @require Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" begin + include("../ext/ManifoldsDistributionsExt/ManifoldsDistributionsExt.jl") + end + + @require OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" begin + @require DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def" begin + include("../ext/ManifoldsOrdinaryDiffEqDiffEqCallbacksExt.jl") + end + end + end + + @require HybridArrays = "1baab800-613f-4b0a-84e4-9cd3431bfbb9" begin + include("../ext/ManifoldsHybridArraysExt.jl") + end end return nothing @@ -673,7 +700,6 @@ export Euclidean, SymmetricPositiveDefinite, SPDFixedDeterminant, SymmetricPositiveSemidefiniteFixedRank, - Symplectic, SymplecticGrassmann, SymplecticMatrices, SymplecticStiefel, @@ -715,11 +741,10 @@ export AbstractPowerManifold, QuotientManifold export ProductManifold, EmbeddedManifold export GraphManifold, GraphManifoldType, VertexManifold, EdgeManifold -export ArrayPartition -export ProjectedPointDistribution, TangentBundle +export TangentBundle export TangentSpace, VectorSpaceFiber, VectorSpaceType, VectorBundle export AbstractVectorTransportMethod, - DifferentiatedRetractionVectorTransport, ParallelTransport, ProjectedPointDistribution + DifferentiatedRetractionVectorTransport, ParallelTransport export PoleLadderTransport, SchildsLadderTransport export ProductVectorTransport export AbstractAffineConnection, ConnectionManifold, LeviCivitaConnection @@ -883,7 +908,6 @@ export ×, minkowski_metric, moment, norm, - normal_tvector_distribution, number_eltype, number_of_coordinates, one, @@ -896,7 +920,6 @@ export ×, parallel_transport_to!, project, project!, - projected_distribution, rand, rand!, real_dimension, @@ -930,7 +953,6 @@ export ×, submanifold, submanifold_component, submanifold_components, - uniform_distribution, var, vector_space_dimension, vector_transport_along, @@ -951,15 +973,18 @@ export ×, # Lie group types & functions export AbstractGroupAction, AbstractGroupOperation, + AbstractGroupVectorRepresentation, ActionDirection, AdditionOperation, CircleGroup, GeneralLinear, GroupManifold, GroupOperationAction, + HybridTangentRepresentation, Identity, LeftAction, LeftInvariantMetric, + LeftInvariantRepresentation, LeftSide, MultiplicationOperation, Orthogonal, @@ -978,6 +1003,7 @@ export AbstractGroupAction, SpecialLinear, SpecialOrthogonal, SpecialUnitary, + TangentVectorRepresentation, TranslationGroup, TranslationAction, Unitary @@ -1005,6 +1031,8 @@ export adjoint_action, compose!, direction, direction_and_side, + exp_inv, + exp_inv!, exp_lie, exp_lie!, group_manifold, @@ -1039,6 +1067,8 @@ export adjoint_action, inverse_translate_diff!, lie_bracket, lie_bracket!, + log_inv, + log_inv!, log_lie, log_lie!, optimal_alignment, diff --git a/src/deprecated.jl b/src/deprecated.jl index eea4fe49ac..8b13789179 100644 --- a/src/deprecated.jl +++ b/src/deprecated.jl @@ -1,9 +1 @@ -@deprecate default_estimation_method(M::AbstractManifold, f) default_approximation_method( - M, - f, -) -@deprecate ExtrinsicEstimation() ExtrinsicEstimation(EfficientEstimator()) - -Base.@deprecate_binding Symplectic SymplecticMatrices -Base.@deprecate_binding SymplecticMatrix SymplecticElement diff --git a/src/distributions.jl b/src/distributions.jl deleted file mode 100644 index c7d361abc9..0000000000 --- a/src/distributions.jl +++ /dev/null @@ -1,67 +0,0 @@ -""" - FVectorvariate - -Structure that subtypes `VariateForm`, indicating that a single sample -is a vector from a fiber of a vector bundle. -""" -struct FVectorvariate <: VariateForm end - -""" - FVectorSupport(space::AbstractManifold, VectorSpaceFiber) - -Value support for vector bundle fiber-valued distributions (values from a fiber of a vector -bundle at a `point` from the given manifold). -For example used for tangent vector-valued distributions. -""" -struct FVectorSupport{TSpace<:VectorSpaceFiber} <: ValueSupport - space::TSpace -end - -""" - FVectorDistribution{TSpace<:VectorSpaceFiber, T} - -An abstract distribution for vector bundle fiber-valued distributions (values from a fiber -of a vector bundle at point `x` from the given manifold). -For example used for tangent vector-valued distributions. -""" -abstract type FVectorDistribution{TSpace<:VectorSpaceFiber} <: - Distribution{FVectorvariate,FVectorSupport{TSpace}} end - -""" - MPointvariate - -Structure that subtypes `VariateForm`, indicating that a single sample -is a point on a manifold. -""" -struct MPointvariate <: VariateForm end - -""" - MPointSupport(M::AbstractManifold) - -Value support for manifold-valued distributions (values from given -[`AbstractManifold`](@extref `ManifoldsBase.AbstractManifold`) `M`). -""" -struct MPointSupport{TM<:AbstractManifold} <: ValueSupport - manifold::TM -end - -""" - MPointDistribution{TM<:AbstractManifold} - -An abstract distribution for points on manifold of type `TM`. -""" -abstract type MPointDistribution{TM<:AbstractManifold} <: - Distribution{MPointvariate,MPointSupport{TM}} end - -""" - support(d::FVectorDistribution) - -Get the object of type `FVectorSupport` for the distribution `d`. -""" -function Distributions.support(::T) where {T<:FVectorDistribution} - return error("support not implemented for type $T") -end - -function uniform_distribution(M::AbstractManifold) - return uniform_distribution(M, allocate_result(M, uniform_distribution)) -end diff --git a/src/groups/GroupManifold.jl b/src/groups/GroupManifold.jl index 28549f33d7..76474a05bd 100644 --- a/src/groups/GroupManifold.jl +++ b/src/groups/GroupManifold.jl @@ -13,22 +13,34 @@ Group manifolds by default forward metric-related operations to the wrapped mani Define the group operation `op` acting on the manifold `manifold`, hence if `op` acts smoothly, this forms a Lie group. """ -struct GroupManifold{𝔽,M<:AbstractManifold{𝔽},O<:AbstractGroupOperation} <: - AbstractDecoratorManifold{𝔽} +struct GroupManifold{ + 𝔽, + M<:AbstractManifold{𝔽}, + O<:AbstractGroupOperation, + VR<:AbstractGroupVectorRepresentation, +} <: AbstractDecoratorManifold{𝔽} manifold::M op::O + vectors::VR end +""" + vector_representation(M::GroupManifold) + +Get the [`AbstractGroupVectorRepresentation`](@ref) of [`GroupManifold`](@ref) `M`. +""" +vector_representation(M::GroupManifold) = M.vectors + @inline function active_traits(f, M::GroupManifold, args...) return merge_traits( - IsGroupManifold(M.op), + IsGroupManifold(M.op, M.vectors), active_traits(f, M.manifold, args...), IsExplicitDecorator(), ) end @inline function active_traits(f, ::AbstractRNG, M::GroupManifold, args...) return merge_traits( - IsGroupManifold(M.op), + IsGroupManifold(M.op, M.vectors), active_traits(f, M.manifold, args...), IsExplicitDecorator(), ) @@ -40,9 +52,17 @@ end decorated_manifold(G::GroupManifold) = G.manifold -(op::AbstractGroupOperation)(M::AbstractManifold) = GroupManifold(M, op) -function (::Type{T})(M::AbstractManifold) where {T<:AbstractGroupOperation} - return GroupManifold(M, T()) +function (op::AbstractGroupOperation)( + M::AbstractManifold, + vectors::AbstractGroupVectorRepresentation, +) + return GroupManifold(M, op, vectors) +end +function (::Type{T})( + M::AbstractManifold, + vectors::AbstractGroupVectorRepresentation, +) where {T<:AbstractGroupOperation} + return GroupManifold(M, T(), vectors) end function inverse_retract( @@ -114,8 +134,8 @@ end @doc raw""" rand(::GroupManifold; vector_at=nothing, σ=1.0) rand!(::GroupManifold, pX; vector_at=nothing, kwargs...) - rand(::TraitList{IsGroupManifold}, M; vector_at=nothing, σ=1.0) - rand!(TraitList{IsGroupManifold}, M, pX; vector_at=nothing, kwargs...) + rand(::TraitList{<:IsGroupManifold}, M; vector_at=nothing, σ=1.0) + rand!(TraitList{<:IsGroupManifold}, M, pX; vector_at=nothing, kwargs...) Compute a random point or tangent vector on a Lie group. diff --git a/src/groups/addition_operation.jl b/src/groups/addition_operation.jl index 885974cf53..57a8190733 100644 --- a/src/groups/addition_operation.jl +++ b/src/groups/addition_operation.jl @@ -22,9 +22,16 @@ Base.:*(e::Identity{AdditionOperation}, ::Identity{AdditionOperation}) = e const AdditionGroupTrait = TraitList{<:IsGroupManifold{AdditionOperation}} -adjoint_action(::AdditionGroupTrait, G::AbstractDecoratorManifold, p, X) = X +adjoint_action(::AdditionGroupTrait, G::AbstractDecoratorManifold, p, X, ::LeftAction) = X -function adjoint_action!(::AdditionGroupTrait, G::AbstractDecoratorManifold, Y, p, X) +function adjoint_action!( + ::AdditionGroupTrait, + G::AbstractDecoratorManifold, + Y, + p, + X, + ::LeftAction, +) return copyto!(G, Y, p, X) end @@ -176,15 +183,3 @@ function translate_diff( ) return X end - -function translate_diff!( - ::AdditionGroupTrait, - G::AbstractDecoratorManifold, - Y, - p, - q, - X, - ::ActionDirectionAndSide, -) - return copyto!(G, Y, p, X) -end diff --git a/src/groups/circle_group.jl b/src/groups/circle_group.jl index 34c97fb430..3bd4d2f3c1 100644 --- a/src/groups/circle_group.jl +++ b/src/groups/circle_group.jl @@ -4,17 +4,27 @@ The circle group is the complex circle ([`Circle(ℂ)`](@ref)) equipped with the group operation of complex multiplication ([`MultiplicationOperation`](@ref)). """ -const CircleGroup = GroupManifold{ℂ,Circle{ℂ},MultiplicationOperation} - -CircleGroup() = GroupManifold(Circle{ℂ}(), MultiplicationOperation()) +const CircleGroup = + GroupManifold{ℂ,Circle{ℂ},MultiplicationOperation,TangentVectorRepresentation} + +function CircleGroup() + return GroupManifold( + Circle{ℂ}(), + MultiplicationOperation(), + TangentVectorRepresentation(), + ) +end @inline function active_traits(f, M::CircleGroup, args...) if is_metric_function(f) #pass to Euclidean by default - but keep Group Decorator for the retraction - return merge_traits(IsGroupManifold(M.op), IsExplicitDecorator()) + return merge_traits( + IsGroupManifold(M.op, TangentVectorRepresentation()), + IsExplicitDecorator(), + ) else return merge_traits( - IsGroupManifold(M.op), + IsGroupManifold(M.op, TangentVectorRepresentation()), IsDefaultMetric(EuclideanMetric()), active_traits(f, M.manifold, args...), IsExplicitDecorator(), #pass to Euclidean by default/last fallback @@ -24,9 +34,13 @@ end Base.show(io::IO, ::CircleGroup) = print(io, "CircleGroup()") -adjoint_action(::CircleGroup, p, X) = X +adjoint_action(::CircleGroup, p, X, ::LeftAction) = X +adjoint_action(::CircleGroup, ::Identity, X, ::LeftAction) = X +adjoint_action(::CircleGroup, p, X, ::RightAction) = X +adjoint_action(::CircleGroup, ::Identity, X, ::RightAction) = X -adjoint_action!(::CircleGroup, Y, p, X) = copyto!(Y, X) +adjoint_action!(::CircleGroup, Y, p, X, ::LeftAction) = copyto!(Y, X) +adjoint_action!(::CircleGroup, Y, p, X, ::RightAction) = copyto!(Y, X) function compose( ::MultiplicationGroupTrait, @@ -77,21 +91,37 @@ lie_bracket(::CircleGroup, X, Y) = zero(X) lie_bracket!(::CircleGroup, Z, X, Y) = fill!(Z, 0) -function translate_diff(::GT, p, q, X, ::ActionDirectionAndSide) where {GT<:CircleGroup} - return map(*, p, X) -end -function translate_diff( - ::CircleGroup, - ::Identity{MultiplicationOperation}, - q, - X, - ::ActionDirectionAndSide, -) - return X +translate_diff(::CircleGroup, p, q, X) = map(*, p, X) +translate_diff(::CircleGroup, p::Identity{MultiplicationOperation}, q, X) = X +for AD in [LeftForwardAction, RightForwardAction, LeftBackwardAction, RightBackwardAction] + @eval begin + function translate_diff(G::CircleGroup, p, q, X, ::$AD) + return translate_diff(G, p, q, X) + end + function translate_diff( + ::CircleGroup, + ::Identity{MultiplicationOperation}, + q, + X, + ::$AD, + ) + return X + end + end end -function translate_diff!(G::CircleGroup, Y, p, q, X, conv::ActionDirectionAndSide) - return copyto!(Y, translate_diff(G, p, q, X, conv)) +_common_translate_diff!(G, Y, p, q, X, conv) = copyto!(Y, translate_diff(G, p, q, X, conv)) +function translate_diff!(G::CircleGroup, Y, p, q, X, conv::LeftForwardAction) + return _common_translate_diff!(G, Y, p, q, X, conv) +end +function translate_diff!(G::CircleGroup, Y, p, q, X, conv::RightForwardAction) + return _common_translate_diff!(G, Y, p, q, X, conv) +end +function translate_diff!(G::CircleGroup, Y, p, q, X, conv::LeftBackwardAction) + return _common_translate_diff!(G, Y, p, q, X, conv) +end +function translate_diff!(G::CircleGroup, Y, p, q, X, conv::RightBackwardAction) + return _common_translate_diff!(G, Y, p, q, X, conv) end function exp_lie(::CircleGroup, X) @@ -124,17 +154,20 @@ end The real circle group is the real circle ([`Circle(ℝ)`](@ref)) equipped with the group operation of addition ([`AdditionOperation`](@ref)). """ -const RealCircleGroup = GroupManifold{ℝ,Circle{ℝ},AdditionOperation} +const RealCircleGroup = + GroupManifold{ℝ,Circle{ℝ},AdditionOperation,LeftInvariantRepresentation} -RealCircleGroup() = GroupManifold(Circle{ℝ}(), AdditionOperation()) +function RealCircleGroup() + return GroupManifold(Circle{ℝ}(), AdditionOperation(), LeftInvariantRepresentation()) +end @inline function active_traits(f, M::RealCircleGroup, args...) if is_metric_function(f) #pass to Euclidean by default - but keep Group Decorator for the retraction - return merge_traits(IsGroupManifold(M.op), IsExplicitDecorator()) + return merge_traits(IsGroupManifold(M.op, M.vectors), IsExplicitDecorator()) else return merge_traits( - IsGroupManifold(M.op), + IsGroupManifold(M.op, M.vectors), HasBiinvariantMetric(), IsDefaultMetric(EuclideanMetric()), active_traits(f, M.manifold, args...), @@ -143,6 +176,17 @@ RealCircleGroup() = GroupManifold(Circle{ℝ}(), AdditionOperation()) end end +adjoint_action(::RealCircleGroup, p, X, ::LeftAction) = X +adjoint_action(::RealCircleGroup, p, X, ::RightAction) = X +adjoint_action(::RealCircleGroup, ::Identity, X, ::LeftAction) = X +adjoint_action(::RealCircleGroup, ::Identity, X, ::RightAction) = X + +for AD in [LeftAction, RightAction] + @eval begin + adjoint_action!(::RealCircleGroup, Y, p, X, ::$AD) = copyto!(Y, X) + end +end + Base.show(io::IO, ::RealCircleGroup) = print(io, "RealCircleGroup()") is_default_metric(::RealCircleGroup, ::EuclideanMetric) = true @@ -198,3 +242,21 @@ function exp_lie(::RealCircleGroup, X) end exp_lie!(::RealCircleGroup, q, X) = (q .= sym_rem(X)) + +translate_diff(::RealCircleGroup, p, q, X) = X +for AD in [LeftForwardAction, RightForwardAction, LeftBackwardAction, RightBackwardAction] + @eval begin + function translate_diff(::RealCircleGroup, p, q, X, ::$AD) + return X + end + function translate_diff( + ::RealCircleGroup, + ::Identity{AdditionOperation}, + q, + X, + ::$AD, + ) + return X + end + end +end diff --git a/src/groups/general_linear.jl b/src/groups/general_linear.jl index 0258f9196b..09f47e8bd7 100644 --- a/src/groups/general_linear.jl +++ b/src/groups/general_linear.jl @@ -21,7 +21,7 @@ end function active_traits(f, ::GeneralLinear, args...) return merge_traits( - IsGroupManifold(MultiplicationOperation()), + IsGroupManifold(MultiplicationOperation(), LeftInvariantRepresentation()), IsEmbeddedManifold(), HasLeftInvariantMetric(), IsDefaultMetric(EuclideanMetric()), @@ -275,9 +275,6 @@ function Base.show(io::IO, M::GeneralLinear{Tuple{Int},𝔽}) where {𝔽} return print(io, "GeneralLinear($n, $(𝔽); parameter=:field)") end -translate_diff(::GeneralLinear, p, q, X, ::LeftForwardAction) = X -translate_diff(::GeneralLinear, p, q, X, ::RightBackwardAction) = p \ X * p - -function translate_diff!(G::GeneralLinear, Y, p, q, X, conv::ActionDirectionAndSide) - return copyto!(Y, translate_diff(G, p, q, X, conv)) -end +# note: this implementation is not optimal +adjoint_action!(::GeneralLinear, Y, p, X, ::LeftAction) = copyto!(Y, p * X * inv(p)) +adjoint_action!(::GeneralLinear, Y, p, X, ::RightAction) = copyto!(Y, p \ X * p) diff --git a/src/groups/general_unitary_groups.jl b/src/groups/general_unitary_groups.jl index b9ebcc0906..c4896767c0 100644 --- a/src/groups/general_unitary_groups.jl +++ b/src/groups/general_unitary_groups.jl @@ -8,16 +8,18 @@ struct GeneralUnitaryMultiplicationGroup{T,𝔽,S} <: AbstractDecoratorManifold{ manifold::GeneralUnitaryMatrices{T,𝔽,S} end +vector_representation(::GeneralUnitaryMultiplicationGroup) = LeftInvariantRepresentation() + @inline function active_traits(f, ::GeneralUnitaryMultiplicationGroup, args...) if is_metric_function(f) #pass to Rotations by default - but keep Group Decorator for the retraction return merge_traits( - IsGroupManifold(MultiplicationOperation()), + IsGroupManifold(MultiplicationOperation(), LeftInvariantRepresentation()), IsExplicitDecorator(), ) else return merge_traits( - IsGroupManifold(MultiplicationOperation()), + IsGroupManifold(MultiplicationOperation(), LeftInvariantRepresentation()), HasBiinvariantMetric(), IsDefaultMetric(EuclideanMetric()), IsExplicitDecorator(), #pass to the inner M by default/last fallback @@ -287,6 +289,15 @@ function manifold_volume(M::GeneralUnitaryMultiplicationGroup) return manifold_volume(M.manifold) end +function Random.rand!(G::GeneralUnitaryMultiplicationGroup, pX; kwargs...) + rand!(G.manifold, pX; kwargs...) + return pX +end +function Random.rand!(rng::AbstractRNG, G::GeneralUnitaryMultiplicationGroup, pX; kwargs...) + rand!(rng, G.manifold, pX; kwargs...) + return pX +end + function translate_diff!( G::GeneralUnitaryMultiplicationGroup, Y, @@ -330,6 +341,14 @@ function translate_diff!( return copyto!(G, Y, inv(G, p) * X * p) end +function adjoint_action!(G::GeneralUnitaryMultiplicationGroup, Y, p, X, ::LeftAction) + copyto!(G, Y, p * X * inv(G, p)) + return Y +end +function adjoint_action!(G::GeneralUnitaryMultiplicationGroup, Y, p, X, ::RightAction) + return copyto!(G, Y, inv(G, p) * X * p) +end + function volume_density(M::GeneralUnitaryMultiplicationGroup, p, X) return volume_density(M.manifold, p, X) end diff --git a/src/groups/group.jl b/src/groups/group.jl index b5e0c59232..537e504362 100644 --- a/src/groups/group.jl +++ b/src/groups/group.jl @@ -19,22 +19,49 @@ For a concrete case the concrete wrapper [`GroupManifold`](@ref) can be used. """ abstract type AbstractGroupOperation end +""" + abstract type AbstractGroupVectorRepresentation end + +An abstract supertype for indicating representation of tangent vectors on a group manifold. +The most common representations are [`LeftInvariantRepresentation`](@ref), +[`TangentVectorRepresentation`](@ref) and [`HybridTangentRepresentation`](@ref). +""" +abstract type AbstractGroupVectorRepresentation end + +""" + TangentVectorRepresentation + +Specify that tangent vectors in a group are stored in a non-invariant way, corresponding to +the storage implied by the underlying manifold. +""" +struct TangentVectorRepresentation <: AbstractGroupVectorRepresentation end + +""" + LeftInvariantRepresentation + +Specify that tangent vectors in a group are stored in Lie algebra using left-invariant +representation. +""" +struct LeftInvariantRepresentation <: AbstractGroupVectorRepresentation end + """ IsGroupManifold{O<:AbstractGroupOperation} <: AbstractTrait A trait to declare an [`AbstractManifold`](@extref `ManifoldsBase.AbstractManifold`) as a manifold with group structure with operation of type `O`. -Using this trait you can turn a manifold that you implement _implictly_ into a Lie group. +Using this trait you can turn a manifold that you implement _implicitly_ into a Lie group. If you wish to decorate an existing manifold with one (or different) [`AbstractGroupAction`](@ref)s, see [`GroupManifold`](@ref). # Constructor - IsGroupManifold(op::AbstractGroupOperation) + IsGroupManifold(op::AbstractGroupOperation, vectors::AbstractGroupVectorRepresentation) """ -struct IsGroupManifold{O<:AbstractGroupOperation} <: AbstractTrait +struct IsGroupManifold{O<:AbstractGroupOperation,VR<:AbstractGroupVectorRepresentation} <: + AbstractTrait op::O + vectors::VR end """ @@ -359,7 +386,7 @@ function isapprox( end @inline function isapprox( - ::TraitList{IsGroupManifold{O}}, + ::TraitList{<:IsGroupManifold{O}}, G::AbstractDecoratorManifold, p::Identity{O}, X, @@ -426,12 +453,12 @@ function is_vector( end @doc raw""" - adjoint_action(G::AbstractDecoratorManifold, p, X) + adjoint_action(G::AbstractDecoratorManifold, p, X, dir=LeftAction()) 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 automorphism ``Ψ_p(q) = pqp⁻¹`` at +If `dir` is `LeftAction()`, it is defined as the differential of the group automorphism ``Ψ_p(q) = pqp⁻¹`` at the identity of `G`. The formula reads @@ -440,28 +467,52 @@ The formula reads ```` where ``e`` is the identity element of `G`. +If `dir` is `RightAction()`, then the formula is +````math +\operatorname{Ad}_p(X) = dΨ_{p^{-1}}(e)[X] +```` + 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::AbstractDecoratorManifold, p, X) -@trait_function adjoint_action(G::AbstractDecoratorManifold, p, Xₑ) -function adjoint_action(::TraitList{<:IsGroupManifold}, G::AbstractDecoratorManifold, p, Xₑ) - Xₚ = translate_diff(G, p, Identity(G), Xₑ, LeftForwardAction()) - Y = inverse_translate_diff(G, p, p, Xₚ, RightBackwardAction()) - return Y +adjoint_action(G::AbstractDecoratorManifold, p, X, dir) +@trait_function adjoint_action(G::AbstractDecoratorManifold, p, Xₑ, dir) +@trait_function adjoint_action!(G::AbstractDecoratorManifold, Y, p, Xₑ, dir) +function adjoint_action( + ::TraitList{<:IsGroupManifold}, + G::AbstractDecoratorManifold, + p, + Xₑ, + dir, +) + BG = base_group(G) + Y = allocate_result(BG, adjoint_action, Xₑ, p) + return adjoint_action!(BG, Y, p, Xₑ, dir) end - -@trait_function adjoint_action!(G::AbstractDecoratorManifold, Y, p, Xₑ) +function adjoint_action(::AbstractDecoratorManifold, ::Identity, Xₑ, ::LeftAction) + return Xₑ +end +function adjoint_action(::AbstractDecoratorManifold, ::Identity, Xₑ, ::RightAction) + return Xₑ +end +# backward compatibility +function adjoint_action(G::AbstractDecoratorManifold, p, X) + return adjoint_action(G, p, X, LeftAction()) +end +function adjoint_action!(G::AbstractDecoratorManifold, Y, p, X) + return adjoint_action!(G, Y, p, X, LeftAction()) +end +# fall back method: the right action is defined from the left action function adjoint_action!( ::TraitList{<:IsGroupManifold}, G::AbstractDecoratorManifold, Y, p, - Xₑ, + X, + ::RightAction, ) - Xₚ = translate_diff(G, p, Identity(G), Xₑ, LeftForwardAction()) - inverse_translate_diff!(G, Y, p, p, Xₚ, RightBackwardAction()) - return Y + BG = base_group(G) + return adjoint_action!(BG, Y, inv(BG, p), X, LeftAction()) end @doc raw""" @@ -525,7 +576,7 @@ function Base.inv(::TraitList{<:IsGroupManifold}, G::AbstractDecoratorManifold, end function Base.inv( - ::TraitList{IsGroupManifold{O}}, + ::TraitList{<:IsGroupManifold{O}}, ::AbstractDecoratorManifold, e::Identity{O}, ) where {O<:AbstractGroupOperation} @@ -535,7 +586,7 @@ end @trait_function inv!(G::AbstractDecoratorManifold, q, p) function inv!( - ::TraitList{IsGroupManifold{O}}, + ::TraitList{<:IsGroupManifold{O}}, G::AbstractDecoratorManifold, q, ::Identity{O}, @@ -544,7 +595,7 @@ function inv!( return identity_element!(BG, q) end function inv!( - ::TraitList{IsGroupManifold{O}}, + ::TraitList{<:IsGroupManifold{O}}, G::AbstractDecoratorManifold, ::Identity{O}, e::Identity{O}, @@ -557,19 +608,29 @@ end Compute the value of differential of inverse ``p^{-1} ∈ \mathcal{G}`` of an element ``p ∈ \mathcal{G}`` at tangent vector `X` at `p`. The result is a tangent vector at ``p^{-1}``. + +*Note*: the default implementation of `inv_diff` and `inv_diff!` +assumes that the tangent vector ``X`` is stored at +the point ``p ∈ \mathcal{G}`` as the vector ``Y ∈ \mathfrak{g}`` + where ``X = pY``. """ inv_diff(G::AbstractDecoratorManifold, p) @trait_function inv_diff(G::AbstractDecoratorManifold, p, X) function inv_diff(::TraitList{<:IsGroupManifold}, G::AbstractDecoratorManifold, p, X) - Y = allocate_result(G, inv_diff, X, p) - return inv_diff!(G, Y, p, X) + return -adjoint_action(base_group(G), p, X) end @trait_function inv_diff!(G::AbstractDecoratorManifold, Y, p, X) +function inv_diff!(::TraitList{<:IsGroupManifold}, G::AbstractDecoratorManifold, Y, p, X) + adjoint_action!(G, Y, p, X) + Y .*= -1 + return Y +end + function Base.copyto!( - ::TraitList{IsGroupManifold{O}}, + ::TraitList{<:IsGroupManifold{O}}, ::AbstractDecoratorManifold, e::Identity{O}, ::Identity{O}, @@ -577,7 +638,7 @@ function Base.copyto!( return e end function Base.copyto!( - ::TraitList{IsGroupManifold{O}}, + ::TraitList{<:IsGroupManifold{O}}, G::AbstractDecoratorManifold, p, ::Identity{O}, @@ -686,7 +747,7 @@ vector to an array representation. The [`vee`](@ref) map is the `hat` map's inverse. """ function hat( - ::TraitList{IsGroupManifold{O}}, + ::TraitList{<:IsGroupManifold{O}}, M::AbstractDecoratorManifold, ::Identity{O}, X, @@ -694,7 +755,7 @@ function hat( return get_vector_lie(M, X, VeeOrthogonalBasis()) end function hat!( - ::TraitList{IsGroupManifold{O}}, + ::TraitList{<:IsGroupManifold{O}}, M::AbstractDecoratorManifold, Y, ::Identity{O}, @@ -728,7 +789,7 @@ vector to a vector representation. The [`hat`](@ref) map is the `vee` map's inverse. """ function vee( - ::TraitList{IsGroupManifold{O}}, + ::TraitList{<:IsGroupManifold{O}}, M::AbstractDecoratorManifold, ::Identity{O}, X, @@ -736,7 +797,7 @@ function vee( return get_coordinates_lie(M, X, VeeOrthogonalBasis()) end function vee!( - ::TraitList{IsGroupManifold{O}}, + ::TraitList{<:IsGroupManifold{O}}, M::AbstractDecoratorManifold, Y, ::Identity{O}, @@ -883,6 +944,12 @@ left or right `conv`ention. The differential transports vectors: ```math (\mathrm{d}τ_p)_q : T_q \mathcal{G} → T_{τ_p q} \mathcal{G}\\ ``` + +*Note*: the default implementation of `translate_diff` and `translate_diff!` +assumes that the tangent vector ``X`` is stored at +the point ``p ∈ \mathcal{G}`` as the vector ``Y ∈ \mathfrak{g}`` + where ``X = pY``. +The implementation at `p = Identity` is independent of the storage choice. """ translate_diff(::AbstractDecoratorManifold, ::Any...) @trait_function translate_diff( @@ -914,6 +981,56 @@ end conv::ActionDirectionAndSide=LeftForwardAction(), ) +function translate_diff!( + ::TraitList{<:IsGroupManifold{<:AbstractGroupOperation,LeftInvariantRepresentation}}, + G::AbstractDecoratorManifold, + Y, + ::Any, + ::Any, + X, + ::LeftForwardAction, +) + return copyto!(G, Y, X) +end +function translate_diff!( + ::TraitList{<:IsGroupManifold{<:AbstractGroupOperation,LeftInvariantRepresentation}}, + G::AbstractDecoratorManifold, + Y, + ::Any, + ::Any, + X, + ::RightForwardAction, +) + return copyto!(G, Y, X) +end +function translate_diff!( + ::TraitList{<:IsGroupManifold{<:AbstractGroupOperation,LeftInvariantRepresentation}}, + G::AbstractDecoratorManifold, + Y, + p, + ::Any, + X, + ::LeftBackwardAction, +) + return adjoint_action!(G, Y, p, X, LeftAction()) +end +function translate_diff!( + ::TraitList{<:IsGroupManifold{<:AbstractGroupOperation,LeftInvariantRepresentation}}, + G::AbstractDecoratorManifold, + Y, + p, + ::Any, + X, + ::RightBackwardAction, +) + return adjoint_action!(G, Y, p, X, RightAction()) +end + +translate_diff(::AbstractDecoratorManifold, ::Identity, q, X, ::LeftForwardAction) = X +translate_diff(::AbstractDecoratorManifold, ::Identity, q, X, ::RightForwardAction) = X +translate_diff(::AbstractDecoratorManifold, ::Identity, q, X, ::LeftBackwardAction) = X +translate_diff(::AbstractDecoratorManifold, ::Identity, q, X, ::RightBackwardAction) = X + @doc raw""" inverse_translate_diff(G::AbstractDecoratorManifold, p, q, X, conv::ActionDirectionAndSide=LeftForwardAction()) @@ -965,6 +1082,44 @@ function inverse_translate_diff!( return translate_diff!(BG, Y, inv(BG, p), q, X, conv) end +""" + log_inv(G::AbstractManifold, p, q) + +Compute logarithmic map on a Lie group `G` invariant to group operation. For groups with a +bi-invariant metric or a Cartan-Schouten connection, this is the same as `log` but for +other groups it may differ. +""" +function log_inv(G::AbstractManifold, p, q) + BG = base_group(G) + return log_lie(BG, compose(BG, inv(BG, p), q)) +end +function log_inv!(G::AbstractManifold, X, p, q) + x = allocate_result(G, inv, p) + BG = base_group(G) + inv!(BG, x, p) + compose!(BG, x, x, q) + log_lie!(BG, X, x) + return X +end + +""" + exp_inv(G::AbstractManifold, p, X, t::Number=1) + +Compute exponential map on a Lie group `G` invariant to group operation. For groups with a +bi-invariant metric or a Cartan-Schouten connection, this is the same as `exp` but for +other groups it may differ. +""" +function exp_inv(G::AbstractManifold, p, X, t::Number=1) + BG = base_group(G) + return compose(BG, p, exp_lie(BG, t * X)) +end +function exp_inv!(G::AbstractManifold, q, p, X) + BG = base_group(G) + exp_lie!(BG, q, X) + compose!(BG, q, p, q) + return q +end + @doc raw""" exp_lie(G, X) exp_lie!(G, q, X) diff --git a/src/groups/heisenberg.jl b/src/groups/heisenberg.jl index 2b9e815def..ddbd1b7a61 100644 --- a/src/groups/heisenberg.jl +++ b/src/groups/heisenberg.jl @@ -26,7 +26,10 @@ function HeisenbergGroup(n::Int; parameter::Symbol=:type) end function active_traits(f, ::HeisenbergGroup, args...) - return merge_traits(IsGroupManifold(MultiplicationOperation()), IsEmbeddedManifold()) + return merge_traits( + IsGroupManifold(MultiplicationOperation(), LeftInvariantRepresentation()), + IsEmbeddedManifold(), + ) end function _heisenberg_a_view(M::HeisenbergGroup, p) @@ -414,8 +417,10 @@ function Base.show(io::IO, M::HeisenbergGroup{Tuple{Int}}) end translate_diff(::HeisenbergGroup, p, q, X, ::LeftForwardAction) = X +translate_diff(::HeisenbergGroup, ::Identity, q, X, ::LeftForwardAction) = X translate_diff(::HeisenbergGroup, p, q, X, ::RightBackwardAction) = p \ X * p +translate_diff(::HeisenbergGroup, ::Identity, q, X, ::RightBackwardAction) = X -function translate_diff!(G::HeisenbergGroup, Y, p, q, X, conv::ActionDirectionAndSide) - return copyto!(Y, translate_diff(G, p, q, X, conv)) -end +# note: this implementation is not optimal +adjoint_action!(::HeisenbergGroup, Y, p, X, ::LeftAction) = copyto!(Y, p * X * inv(p)) +adjoint_action!(::HeisenbergGroup, Y, p, X, ::RightAction) = copyto!(Y, p \ X * p) diff --git a/src/groups/power_group.jl b/src/groups/power_group.jl index 1857963d31..7b9f23d322 100644 --- a/src/groups/power_group.jl +++ b/src/groups/power_group.jl @@ -37,7 +37,7 @@ function PowerGroup(manifold::AbstractPowerManifold) error("All powered manifold must be or decorate a group.") end op = ProductOperation() - return GroupManifold(manifold, op) + return GroupManifold(manifold, op, vector_representation(manifold.manifold)) end function ManifoldsBase._access_nested( @@ -77,10 +77,10 @@ end @inline function active_traits(f, M::PowerGroup, args...) if is_metric_function(f) #pass to manifold by default - but keep Group Decorator for the retraction - return merge_traits(IsGroupManifold(M.op), IsExplicitDecorator()) + return merge_traits(IsGroupManifold(M.op, M.vectors), IsExplicitDecorator()) else return merge_traits( - IsGroupManifold(M.op), + IsGroupManifold(M.op, M.vectors), active_traits(f, M.manifold, args...), IsExplicitDecorator(), ) @@ -259,7 +259,66 @@ function inverse_translate!( return x end -function translate_diff!(G::PowerGroup, Y, p, q, X, conv::ActionDirectionAndSide) +function _common_power_adjoint_action!(G, Y, p, X, conv) + GM = G.manifold + N = GM.manifold + rep_size = representation_size(N) + for i in get_iterator(GM) + adjoint_action!( + N, + _write(GM, rep_size, Y, i), + _read(GM, rep_size, p, i), + _read(GM, rep_size, X, i), + conv, + ) + end + return Y +end +function adjoint_action!(G::PowerGroup, Y, p, X, conv::LeftAction) + return _common_power_adjoint_action!(G, Y, p, X, conv) +end +function adjoint_action!(G::PowerGroup, Y, p, X, conv::RightAction) + return _common_power_adjoint_action!(G, Y, p, X, conv) +end + +function _common_power_replacing_adjoint_action!(G, Y, p, X, conv) + GM = G.manifold + N = GM.manifold + rep_size = representation_size(N) + for i in get_iterator(GM) + Y[i...] = + adjoint_action(N, _read(GM, rep_size, p, i), _read(GM, rep_size, X, i), conv) + end + return Y +end +function adjoint_action!(G::PowerGroupNestedReplacing, Y, p, X, conv::LeftAction) + return _common_power_replacing_adjoint_action!(G, Y, p, X, conv) +end +function adjoint_action!(G::PowerGroupNestedReplacing, Y, p, X, conv::RightAction) + return _common_power_replacing_adjoint_action!(G, Y, p, X, conv) +end + +function translate_diff!(G::PowerGroup, Y, p, q, X, conv::LeftForwardAction) + return _common_power_translate_diff!(G, Y, p, q, X, conv) +end +function translate_diff!(G::PowerGroup, Y, p, q, X, conv::RightForwardAction) + return _common_power_translate_diff!(G, Y, p, q, X, conv) +end +function translate_diff!(G::PowerGroup, Y, p, q, X, conv::LeftBackwardAction) + return _common_power_translate_diff!(G, Y, p, q, X, conv) +end +function translate_diff!(G::PowerGroup, Y, p, q, X, conv::RightBackwardAction) + return _common_power_translate_diff!(G, Y, p, q, X, conv) +end + +function _common_power_translate_diff!( + G::PowerGroup, + Y, + p, + q, + X, + conv::ActionDirectionAndSide, +) GM = G.manifold N = GM.manifold rep_size = representation_size(N) @@ -275,7 +334,27 @@ function translate_diff!(G::PowerGroup, Y, p, q, X, conv::ActionDirectionAndSide end return Y end + +function translate_diff!(G::PowerGroupNestedReplacing, Y, p, q, X, conv::LeftForwardAction) + return _common_power_replacing_translate_diff!(G, Y, p, q, X, conv) +end +function translate_diff!(G::PowerGroupNestedReplacing, Y, p, q, X, conv::RightForwardAction) + return _common_power_replacing_translate_diff!(G, Y, p, q, X, conv) +end +function translate_diff!(G::PowerGroupNestedReplacing, Y, p, q, X, conv::LeftBackwardAction) + return _common_power_replacing_translate_diff!(G, Y, p, q, X, conv) +end function translate_diff!( + G::PowerGroupNestedReplacing, + Y, + p, + q, + X, + conv::RightBackwardAction, +) + return _common_power_replacing_translate_diff!(G, Y, p, q, X, conv) +end +function _common_power_replacing_translate_diff!( G::PowerGroupNestedReplacing, Y, p, diff --git a/src/groups/product_group.jl b/src/groups/product_group.jl index b630edf4ef..ae8bce3a10 100644 --- a/src/groups/product_group.jl +++ b/src/groups/product_group.jl @@ -19,43 +19,59 @@ one. This type is mostly useful for equipping the direct product of group manifo # Constructor ProductGroup(manifold::ProductManifold) """ -function ProductGroup(manifold::ProductManifold{𝔽}) where {𝔽} +function ProductGroup(manifold::ProductManifold, vectors::AbstractGroupVectorRepresentation) if !all(is_group_manifold, manifold.manifolds) error("All submanifolds of product manifold must be or decorate groups.") end op = ProductOperation() - return GroupManifold(manifold, op) + return GroupManifold(manifold, op, vectors) end @inline function active_traits(f, M::ProductGroup, args...) if is_metric_function(f) #pass to manifold by default - but keep Group Decorator for the retraction - return merge_traits(IsGroupManifold(M.op), IsExplicitDecorator()) + return merge_traits(IsGroupManifold(M.op, M.vectors), IsExplicitDecorator()) else return merge_traits( - IsGroupManifold(M.op), + IsGroupManifold(M.op, M.vectors), active_traits(f, M.manifold, args...), IsExplicitDecorator(), ) end end -function adjoint_inv_diff!(G::ProductGroup, Y, p, X) +function _common_product_adjoint_action!(G, Y, p, X, conv) M = G.manifold map( - adjoint_inv_diff!, + adjoint_action!, M.manifolds, submanifold_components(G, Y), submanifold_components(G, p), submanifold_components(G, X), + repeated(conv), ) return Y end -function identity_element(G::ProductGroup) +function adjoint_action!(G::ProductGroup, Y, p, X, conv::LeftAction) + return _common_product_adjoint_action!(G, Y, p, X, conv) +end +function adjoint_action!(G::ProductGroup, Y, p, X, conv::RightAction) + return _common_product_adjoint_action!(G, Y, p, X, conv) +end + +function adjoint_inv_diff!(G::ProductGroup, Y, p, X) M = G.manifold - return ArrayPartition(map(identity_element, M.manifolds)) + map( + adjoint_inv_diff!, + M.manifolds, + submanifold_components(G, Y), + submanifold_components(G, p), + submanifold_components(G, X), + ) + return Y end + function identity_element!(G::ProductGroup, p) pes = submanifold_components(G, p) M = G.manifold @@ -127,16 +143,6 @@ function inv_diff!(G::ProductGroup, Y, p, X) end _compose(G::ProductGroup, p, q) = _compose(G.manifold, p, q) -function _compose(M::ProductManifold, p::ArrayPartition, q::ArrayPartition) - return ArrayPartition( - map( - compose, - M.manifolds, - submanifold_components(M, p), - submanifold_components(M, q), - )..., - ) -end _compose!(G::ProductGroup, x, p, q) = _compose!(G.manifold, x, p, q) function _compose!(M::ProductManifold, x, p, q) @@ -150,23 +156,6 @@ function _compose!(M::ProductManifold, x, p, q) return x end -function translate( - M::ProductGroup, - p::ArrayPartition, - q::ArrayPartition, - conv::ActionDirectionAndSide, -) - return ArrayPartition( - map( - translate, - M.manifold.manifolds, - submanifold_components(M, p), - submanifold_components(M, q), - repeated(conv), - )..., - ) -end - function translate!(M::ProductGroup, x, p, q, conv::ActionDirectionAndSide) map( translate!, @@ -179,19 +168,6 @@ function translate!(M::ProductGroup, x, p, q, conv::ActionDirectionAndSide) return x end -function inverse_translate(G::ProductGroup, p, q, conv::ActionDirectionAndSide) - M = G.manifold - return ArrayPartition( - map( - inverse_translate, - M.manifolds, - submanifold_components(G, p), - submanifold_components(G, q), - repeated(conv), - )..., - ) -end - function inverse_translate!(G::ProductGroup, x, p, q, conv::ActionDirectionAndSide) M = G.manifold map( @@ -205,21 +181,34 @@ function inverse_translate!(G::ProductGroup, x, p, q, conv::ActionDirectionAndSi return x end -function translate_diff(G::ProductGroup, p, q, X, conv::ActionDirectionAndSide) - M = G.manifold - return ArrayPartition( - map( - translate_diff, - M.manifolds, - submanifold_components(G, p), - submanifold_components(G, q), - submanifold_components(G, X), - repeated(conv), - )..., - ) +function _common_product_translate_diff end + +function translate_diff(G::ProductGroup, p, q, X, conv::LeftForwardAction) + return _common_product_translate_diff(G, p, q, X, conv) +end +function translate_diff(G::ProductGroup, p, q, X, conv::RightForwardAction) + return _common_product_translate_diff(G, p, q, X, conv) +end +function translate_diff(G::ProductGroup, p, q, X, conv::LeftBackwardAction) + return _common_product_translate_diff(G, p, q, X, conv) end +function translate_diff(G::ProductGroup, p, q, X, conv::RightBackwardAction) + return _common_product_translate_diff(G, p, q, X, conv) +end + +translate_diff(::ProductGroup, ::Identity, q, X, ::LeftForwardAction) = X +translate_diff(::ProductGroup, ::Identity, q, X, ::RightForwardAction) = X +translate_diff(::ProductGroup, ::Identity, q, X, ::LeftBackwardAction) = X +translate_diff(::ProductGroup, ::Identity, q, X, ::RightBackwardAction) = X -function translate_diff!(G::ProductGroup, Y, p, q, X, conv::ActionDirectionAndSide) +function _common_product_translate_diff!( + G::ProductGroup, + Y, + p, + q, + X, + conv::ActionDirectionAndSide, +) M = G.manifold map( translate_diff!, @@ -233,18 +222,17 @@ function translate_diff!(G::ProductGroup, Y, p, q, X, conv::ActionDirectionAndSi return Y end -function inverse_translate_diff(G::ProductGroup, p, q, X, conv::ActionDirectionAndSide) - M = G.manifold - return ArrayPartition( - map( - inverse_translate_diff, - M.manifolds, - submanifold_components(G, p), - submanifold_components(G, q), - submanifold_components(G, X), - repeated(conv), - )..., - ) +function translate_diff!(G::ProductGroup, Y, p, q, X, conv::LeftForwardAction) + return _common_product_translate_diff!(G, Y, p, q, X, conv) +end +function translate_diff!(G::ProductGroup, Y, p, q, X, conv::RightForwardAction) + return _common_product_translate_diff!(G, Y, p, q, X, conv) +end +function translate_diff!(G::ProductGroup, Y, p, q, X, conv::LeftBackwardAction) + return _common_product_translate_diff!(G, Y, p, q, X, conv) +end +function translate_diff!(G::ProductGroup, Y, p, q, X, conv::RightBackwardAction) + return _common_product_translate_diff!(G, Y, p, q, X, conv) end function inverse_translate_diff!(G::ProductGroup, Y, p, q, X, conv::ActionDirectionAndSide) @@ -261,17 +249,6 @@ function inverse_translate_diff!(G::ProductGroup, Y, p, q, X, conv::ActionDirect return Y end -function Base.exp(M::ProductGroup, p::Identity{ProductOperation}, X::ArrayPartition) - return ArrayPartition( - map( - exp, - M.manifold.manifolds, - submanifold_components(M, p), - submanifold_components(M, X), - )..., - ) -end - function exp!(M::ProductGroup, q, p::Identity{ProductOperation}, X) map( exp!, @@ -283,28 +260,12 @@ function exp!(M::ProductGroup, q, p::Identity{ProductOperation}, X) return q end -function exp_lie(G::ProductGroup, X) - M = G.manifold - return ArrayPartition(map(exp_lie, M.manifolds, submanifold_components(G, X))...) -end - function exp_lie!(G::ProductGroup, q, X) M = G.manifold map(exp_lie!, M.manifolds, submanifold_components(G, q), submanifold_components(G, X)) return q end -function Base.log(M::ProductGroup, p::Identity{ProductOperation}, q::ArrayPartition) - return ArrayPartition( - map( - log, - M.manifold.manifolds, - submanifold_components(M, p), - submanifold_components(M, q), - )..., - ) -end - function log!(M::ProductGroup, X, p::Identity{ProductOperation}, q) map( log!, @@ -330,36 +291,7 @@ function log_lie!(G::ProductGroup, X, q::Identity{ProductOperation}) return X end -Base.@propagate_inbounds function Base.getindex( - p::ArrayPartition, - M::ProductGroup, - i::Union{Integer,Colon,AbstractVector,Val}, -) - return getindex(p, base_manifold(M), i) -end - -Base.@propagate_inbounds function Base.setindex!( - q::ArrayPartition, - p, - M::ProductGroup, - i::Union{Integer,Colon,AbstractVector,Val}, -) - return setindex!(q, p, base_manifold(M), i) -end - -# these isapprox methods are here just to reduce time-to-first-isapprox -function isapprox(G::ProductGroup, p::ArrayPartition, q::ArrayPartition; kwargs...) - return isapprox(G.manifold, p, q; kwargs...) -end -function isapprox( - G::ProductGroup, - p::ArrayPartition, - X::ArrayPartition, - Y::ArrayPartition; - kwargs..., -) - return isapprox(G.manifold, p, X, Y; kwargs...) -end +# this isapprox method is here just to reduce time-to-first-isapprox function isapprox(G::ProductGroup, ::Identity{ProductOperation}, X, Y; kwargs...) return isapprox(G.manifold, identity_element(G), X, Y; kwargs...) end diff --git a/src/groups/rotation_translation_action.jl b/src/groups/rotation_translation_action.jl index 1883198be6..6f881b11fd 100644 --- a/src/groups/rotation_translation_action.jl +++ b/src/groups/rotation_translation_action.jl @@ -39,11 +39,11 @@ Alias for [`RotationTranslationAction`](@ref) where the manifold `M` is [`Euclid or [`TranslationGroup`](@ref) with size of type `TE`, and [`SpecialEuclidean`](@ref) group has size type `TSE`. """ -const RotationTranslationActionOnVector{TAD,𝔽,TE,TSE} = RotationTranslationAction{ +const RotationTranslationActionOnVector{TAD,𝔽,TE,TSE,SE_GVR} = RotationTranslationAction{ TAD, <:Union{Euclidean{TE,𝔽},TranslationGroup{TE,𝔽}}, - SpecialEuclidean{TSE}, -} where {TAD<:ActionDirection,𝔽,TE,TSE} + SpecialEuclidean{TSE,SE_GVR}, +} where {TAD<:ActionDirection,𝔽,TE,TSE,SE_GVR<:AbstractGroupVectorRepresentation} base_group(A::RotationTranslationAction) = A.SEn @@ -53,29 +53,6 @@ function switch_direction(A::RotationTranslationAction{TAD}) where {TAD<:ActionD return RotationTranslationAction(A.manifold, A.SEn, switch_direction(TAD())) end -""" - apply(::RotationTranslationActionOnVector{LeftAction}, a::ArrayPartition, p) - -Rotate point `p` by `a.x[2]` and translate it by `a.x[1]`. -""" -function apply(::RotationTranslationActionOnVector{LeftAction}, a::ArrayPartition, p) - return a.x[2] * p + a.x[1] -end -function apply( - ::RotationTranslationActionOnVector{LeftAction}, - a::SpecialEuclideanIdentity, - p, -) - return p -end -""" - apply(::RotationTranslationActionOnVector{RightAction}, a::ArrayPartition, p) - -Translate point `p` by `-a.x[1]` and rotate it by inverse of `a.x[2]`. -""" -function apply(::RotationTranslationActionOnVector{RightAction}, a::ArrayPartition, p) - return a.x[2] \ (p - a.x[1]) -end function apply( ::RotationTranslationActionOnVector{RightAction}, a::SpecialEuclideanIdentity, @@ -84,11 +61,6 @@ function apply( return p end -function apply!(::RotationTranslationActionOnVector{LeftAction}, q, a::ArrayPartition, p) - mul!(q, a.x[2], p) - q .+= a.x[1] - return q -end function apply!( ::RotationTranslationActionOnVector{LeftAction}, q, @@ -99,50 +71,6 @@ function apply!( return q end -""" - inverse_apply(::RotationTranslationActionOnVector{LeftAction}, a::ArrayPartition, p) - -Translate point `p` by `-a.x[1]` and rotate it by inverse of `a.x[2]`. -""" -function inverse_apply( - ::RotationTranslationActionOnVector{LeftAction}, - a::ArrayPartition, - p, -) - return a.x[2] \ (p - a.x[1]) -end -""" - inverse_apply(::RotationTranslationActionOnVector{RightAction}, a::ArrayPartition, p) - -Rotate point `p` by `a.x[2]` and translate it by `a.x[1]`. -""" -function inverse_apply( - ::RotationTranslationActionOnVector{RightAction}, - a::ArrayPartition, - p, -) - return a.x[2] * p + a.x[1] -end - -""" - apply_diff( - ::RotationTranslationActionOnVector{LeftAction}, - a::ArrayPartition, - p, - X, - ) - -Compute differential of `apply` on left [`RotationTranslationActionOnVector`](@ref), -with respect to `p`, i.e. left-multiply vector `X` tangent at `p` by `a.x[2]`. -""" -function apply_diff( - ::RotationTranslationActionOnVector{LeftAction}, - a::ArrayPartition, - p, - X, -) - return a.x[2] * X -end function apply_diff( ::RotationTranslationActionOnVector{LeftAction}, ::SpecialEuclideanIdentity, @@ -151,25 +79,6 @@ function apply_diff( ) return X end -""" - apply_diff( - ::RotationTranslationActionOnVector{RightAction}, - a::ArrayPartition, - p, - X, - ) - -Compute differential of `apply` on right [`RotationTranslationActionOnVector`](@ref), -with respect to `p`, i.e. left-divide vector `X` tangent at `p` by `a.x[2]`. -""" -function apply_diff( - ::RotationTranslationActionOnVector{RightAction}, - a::ArrayPartition, - p, - X, -) - return a.x[2] \ X -end function apply_diff( ::RotationTranslationActionOnVector{RightAction}, a::SpecialEuclideanIdentity, @@ -179,16 +88,6 @@ function apply_diff( return X end -function apply_diff!( - ::RotationTranslationActionOnVector{LeftAction}, - Y, - a::ArrayPartition, - p, - X, -) - mul!(Y, a.x[2], X) - return Y -end function apply_diff!( ::RotationTranslationActionOnVector{LeftAction}, Y, @@ -198,16 +97,6 @@ function apply_diff!( ) return copyto!(Y, X) end -function apply_diff!( - ::RotationTranslationActionOnVector{RightAction}, - Y, - a::ArrayPartition, - p, - X, -) - Y .= a.x[2] \ X - return Y -end function apply_diff!( ::RotationTranslationActionOnVector{RightAction}, Y, @@ -238,34 +127,6 @@ function apply_diff_group( return X.x[2] * p end -function apply_diff_group!( - ::RotationTranslationActionOnVector{LeftAction}, - Y, - ::SpecialEuclideanIdentity, - X::ArrayPartition, - p, -) - Y .= X.x[2] * p - return Y -end - -function inverse_apply_diff( - ::RotationTranslationActionOnVector{LeftAction}, - a::ArrayPartition, - p, - X, -) - return a.x[2] \ X -end -function inverse_apply_diff( - ::RotationTranslationActionOnVector{RightAction}, - a::ArrayPartition, - p, - X, -) - return a.x[2] * X -end - ### @doc raw""" @@ -306,18 +167,10 @@ end const LeftColumnwiseSpecialEuclideanAction{TM<:AbstractManifold,TSE<:SpecialEuclidean} = ColumnwiseSpecialEuclideanAction{LeftAction,TM,TSE} -function apply(::LeftColumnwiseSpecialEuclideanAction, a::ArrayPartition, p) - return a.x[2] * p .+ a.x[1] -end function apply(::LeftColumnwiseSpecialEuclideanAction, ::SpecialEuclideanIdentity, p) return p end -function apply!(::LeftColumnwiseSpecialEuclideanAction, q, a::ArrayPartition, p) - map((qrow, prow) -> mul!(qrow, a.x[2], prow), eachcol(q), eachcol(p)) - q .+= a.x[1] - return q -end function apply!(::LeftColumnwiseSpecialEuclideanAction, q, a::SpecialEuclideanIdentity, p) copyto!(q, p) return q @@ -326,30 +179,3 @@ end base_group(A::LeftColumnwiseSpecialEuclideanAction) = A.SE group_manifold(A::LeftColumnwiseSpecialEuclideanAction) = A.manifold - -function inverse_apply(::LeftColumnwiseSpecialEuclideanAction, a::ArrayPartition, p) - return a.x[2] \ (p .- a.x[1]) -end - -@doc raw""" - optimal_alignment(A::LeftColumnwiseSpecialEuclideanAction, p, q) - -Compute optimal alignment of `p` to `q` under the forward left [`ColumnwiseSpecialEuclideanAction`](@ref). -The algorithm, in sequence, computes optimal translation and optimal rotation. -""" -function optimal_alignment( - A::LeftColumnwiseSpecialEuclideanAction{<:AbstractManifold,<:SpecialEuclidean}, - p, - q, -) - N = _get_parameter(A.SE) - tr_opt = mean(q; dims=1) - mean(p; dims=1) - p_moved = p .+ tr_opt - - Ostar = optimal_alignment( - ColumnwiseMultiplicationAction(A.manifold, SpecialOrthogonal(N)), - p_moved, - q, - ) - return ArrayPartition(tr_opt, Ostar) -end diff --git a/src/groups/semidirect_product_group.jl b/src/groups/semidirect_product_group.jl index 5d850fed4e..08857be743 100644 --- a/src/groups/semidirect_product_group.jl +++ b/src/groups/semidirect_product_group.jl @@ -13,8 +13,25 @@ function Base.show(io::IO, op::SemidirectProductOperation) return print(io, "SemidirectProductOperation($(op.action))") end -const SemidirectProductGroup{𝔽,N,H,A<:AbstractGroupAction} = - GroupManifold{𝔽,ProductManifold{𝔽,Tuple{N,H}},SemidirectProductOperation{A}} +""" + struct HybridTangentRepresentation <: AbstractGroupVectorRepresentation end + +Tangent vector representation on [`SemidirectProductGroup`](@ref) such as +[`SpecialEuclidean`](@ref) that corresponds to simple product structure of underlying +groups. +""" +struct HybridTangentRepresentation <: AbstractGroupVectorRepresentation end + +const SemidirectProductGroup{ + 𝔽, + N, + H, + A<:AbstractGroupAction, + GVR<:AbstractGroupVectorRepresentation, +} = GroupManifold{𝔽,ProductManifold{𝔽,Tuple{N,H}},SemidirectProductOperation{A},GVR} + +const SemidirectProductGroupHVR{𝔽,N,H,A<:AbstractGroupAction} = + SemidirectProductGroup{𝔽,N,H,A,HybridTangentRepresentation} @doc raw""" SemidirectProductGroup(N::GroupManifold, H::GroupManifold, A::AbstractGroupAction) @@ -38,20 +55,13 @@ function SemidirectProductGroup( N::AbstractDecoratorManifold{𝔽}, H::AbstractDecoratorManifold{𝔽}, A::AbstractGroupAction, + vectors::AbstractGroupVectorRepresentation, ) where {𝔽} N === group_manifold(A) || error("Subgroup $(N) must be the G-manifold of action $(A)") H === base_group(A) || error("Subgroup $(H) must be the base group of action $(A)") op = SemidirectProductOperation(A) M = ProductManifold(N, H) - return GroupManifold(M, op) -end - -function allocate_result(G::SemidirectProductGroup, ::typeof(identity_element)) - M = base_manifold(G) - N, H = M.manifolds - np = allocate_result(N, identity_element) - hp = allocate_result(H, identity_element) - return ArrayPartition(np, hp) + return GroupManifold(M, op, vectors) end """ @@ -136,9 +146,10 @@ function _compose!(G::SemidirectProductGroup, x, p, q) end @doc raw""" - translate_diff(G::SemidirectProductGroup, p, q, X, conX::LeftForwardAction) + translate_diff(G::SemidirectProductGroupHVR, p, q, X, conX::LeftForwardAction) -Perform differential of the left translation on the semidirect product group `G`. +Perform differential of the left translation on the semidirect product group `G` +with `HybridTangentRepresentation`. Since the left translation is defined as (cf. [`SemidirectProductGroup`](@ref)): @@ -152,9 +163,9 @@ then its differential can be computed as \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::LeftForwardAction) +translate_diff(G::SemidirectProductGroupHVR, p, q, X, conX::LeftForwardAction) -function translate_diff!(G::SemidirectProductGroup, Y, p, q, X, conX::LeftForwardAction) +function translate_diff!(G::SemidirectProductGroupHVR, Y, p, q, X, conX::LeftForwardAction) M = base_manifold(G) N, H = M.manifolds A = G.op.action @@ -248,23 +259,6 @@ function isapprox( return isapprox(G, identity_element(G), X, Y; kwargs...) end -Base.@propagate_inbounds function Base.getindex( - p::ArrayPartition, - M::SemidirectProductGroup, - i::Union{Integer,Colon,AbstractVector,Val}, -) - return getindex(p, base_manifold(M), i) -end - -Base.@propagate_inbounds function Base.setindex!( - q::ArrayPartition, - p, - M::SemidirectProductGroup, - i::Union{Integer,Colon,AbstractVector,Val}, -) - return setindex!(q, p, base_manifold(M), i) -end - function submanifold_components( M::ProductManifold, ::Identity{<:SemidirectProductOperation}, diff --git a/src/groups/special_euclidean.jl b/src/groups/special_euclidean.jl index 9382b2236d..d84559948d 100644 --- a/src/groups/special_euclidean.jl +++ b/src/groups/special_euclidean.jl @@ -1,5 +1,8 @@ @doc raw""" - SpecialEuclidean(n) + SpecialEuclidean( + n::Int; + vectors::AbstractGroupVectorRepresentation=LeftInvariantRepresentation() + ) Special Euclidean group ``\mathrm{SE}(n)``, the group of rigid motions. @@ -17,13 +20,20 @@ This constructor is equivalent to calling ```julia Tn = TranslationGroup(n) SOn = SpecialOrthogonal(n) -SemidirectProductGroup(Tn, SOn, RotationAction(Tn, SOn)) +SemidirectProductGroup(Tn, SOn, RotationAction(Tn, SOn), vectors) ``` Points on ``\mathrm{SE}(n)`` may be represented as points on the underlying product manifold ``\mathrm{T}(n) × \mathrm{SO}(n)``. For group-specific functions, they may also be represented as affine matrices with size `(n + 1, n + 1)` (see [`affine_matrix`](@ref)), for which the group operation is [`MultiplicationOperation`](@ref). + +There are two supported conventions for tangent vector storage, which can be selected +using the `vectors` keyword argument: +* [`LeftInvariantRepresentation`](@ref) (default one), which corresponds to left-invariant + storage commonly used in other Lie groups. +* [`HybridTangentRepresentation`](@ref) which corresponds to the representation implied by + product manifold structure of underlying groups. """ const SpecialEuclidean{T} = SemidirectProductGroup{ ℝ, @@ -35,11 +45,15 @@ const SpecialEuclidean{T} = SemidirectProductGroup{ const SpecialEuclideanManifold{N} = ProductManifold{ℝ,Tuple{TranslationGroup{N,ℝ},SpecialOrthogonal{N}}} -function SpecialEuclidean(n; parameter::Symbol=:type) +function SpecialEuclidean( + n::Int; + vectors::AbstractGroupVectorRepresentation=LeftInvariantRepresentation(), + parameter::Symbol=:type, +) Tn = TranslationGroup(n; parameter=parameter) SOn = SpecialOrthogonal(n; parameter=parameter) A = RotationAction(Tn, SOn) - return SemidirectProductGroup(Tn, SOn, A) + return SemidirectProductGroup(Tn, SOn, A, vectors) end const SpecialEuclideanOperation{N} = SemidirectProductOperation{ @@ -47,16 +61,24 @@ const SpecialEuclideanOperation{N} = SemidirectProductOperation{ } const SpecialEuclideanIdentity{N} = Identity{SpecialEuclideanOperation{N}} -function Base.show(io::IO, ::SpecialEuclidean{TypeParameter{Tuple{n}}}) where {n} - return print(io, "SpecialEuclidean($(n))") +function Base.show(io::IO, G::SpecialEuclidean{TypeParameter{Tuple{n}}}) where {n} + if vector_representation(G) isa LeftInvariantRepresentation + return print(io, "SpecialEuclidean($(n))") + else + return print(io, "SpecialEuclidean($(n); vectors=$(G.vectors))") + end end function Base.show(io::IO, G::SpecialEuclidean{Tuple{Int}}) n = _get_parameter(G) - return print(io, "SpecialEuclidean($(n); parameter=:field)") + if vector_representation(G) isa LeftInvariantRepresentation + return print(io, "SpecialEuclidean($(n); parameter=:field)") + else + return print(io, "SpecialEuclidean($(n); parameter=:field, vectors=$(G.vectors))") + end end @inline function active_traits(f, M::SpecialEuclidean, args...) - return merge_traits(IsGroupManifold(M.op), IsExplicitDecorator()) + return merge_traits(IsGroupManifold(M.op, M.vectors), IsExplicitDecorator()) end """ @@ -142,7 +164,11 @@ Base.@propagate_inbounds function _padvector!( end @doc raw""" - adjoint_action(::SpecialEuclidean{TypeParameter{Tuple{3}}}, p, fX::TFVector{<:Any,VeeOrthogonalBasis{ℝ}}) + adjoint_action( + ::SpecialEuclidean{TypeParameter{Tuple{3}},<:HybridTangentRepresentation}, + p, + fX::TFVector{<:Any,VeeOrthogonalBasis{ℝ}}, + ) Adjoint action of the [`SpecialEuclidean`](@ref) group on the vector with coefficients `fX` tangent at point `p`. @@ -153,7 +179,7 @@ matrix part of `p`, `r` is the translation part of `fX` and `ω` is the rotation ``×`` is the cross product and ``⋅`` is the matrix product. """ function adjoint_action( - ::SpecialEuclidean{TypeParameter{Tuple{3}}}, + ::SpecialEuclidean{TypeParameter{Tuple{3}},<:HybridTangentRepresentation}, p, fX::TFVector{<:Any,VeeOrthogonalBasis{ℝ}}, ) @@ -600,11 +626,6 @@ algebra. For the matrix representation (which can be obtained using [`screw_matr the formula is ``[X, Y] = XY-YX``, while in the `ArrayPartition` 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::ArrayPartition, Y::ArrayPartition) - nX, hX = submanifold_components(G, X) - nY, hY = submanifold_components(G, Y) - return ArrayPartition(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 @@ -633,7 +654,14 @@ is the translation part of `p` and ``X_t`` is the translation part of `X`. """ translate_diff(G::SpecialEuclidean, p, q, X, ::RightBackwardAction) -function translate_diff!(G::SpecialEuclidean, Y, p, q, X, ::RightBackwardAction) +function translate_diff!( + G::SpecialEuclidean{T,<:HybridTangentRepresentation}, + Y, + p, + q, + X, + ::RightBackwardAction, +) where {T} np, hp = submanifold_components(G, p) nq, hq = submanifold_components(G, q) nX, hX = submanifold_components(G, X) @@ -644,6 +672,19 @@ function translate_diff!(G::SpecialEuclidean, Y, p, q, X, ::RightBackwardAction) return Y end +function adjoint_action!(G::SpecialEuclidean, Y, p, Xₑ, ::LeftAction) + np, hp = submanifold_components(G, p) + n, h = submanifold_components(G, Y) + nX, hX = submanifold_components(G, Xₑ) + H = submanifold(G, 2) + adjoint_action!(H, h, hp, hX, LeftAction()) + A = G.op.action + apply!(A, n, hp, nX) + LinearAlgebra.axpy!(-1, apply_diff_group(A, Identity(H), h, np), n) + @inbounds _padvector!(G, Y) + return Y +end + @doc raw""" SpecialEuclideanInGeneralLinear @@ -653,13 +694,22 @@ Note that this is *not* a transparently isometric embedding. # Constructor - SpecialEuclideanInGeneralLinear(n) + SpecialEuclideanInGeneralLinear( + n::Int; + se_vectors::AbstractGroupVectorRepresentation=LeftInvariantVectorRepresentation(), + ) + +Where `se_vectors` is the tangent vector representation of the [`SpecialEuclidean`](@ref) +group to be used. """ const SpecialEuclideanInGeneralLinear = EmbeddedManifold{ℝ,<:SpecialEuclidean,<:GeneralLinear} -function SpecialEuclideanInGeneralLinear(n) - return EmbeddedManifold(SpecialEuclidean(n), GeneralLinear(n + 1)) +function SpecialEuclideanInGeneralLinear( + n::Int; + se_vectors::AbstractGroupVectorRepresentation=LeftInvariantVectorRepresentation(), +) + return EmbeddedManifold(SpecialEuclidean(n; vectors=se_vectors), GeneralLinear(n + 1)) end """ @@ -675,10 +725,9 @@ end """ embed(M::SpecialEuclideanInGeneralLinear, p, X) -Embed the tangent vector X at point `p` on [`SpecialEuclidean`](@ref) in the +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. +`SpecialEuclidean`. The embedding is similar from the one defined by [`screw_matrix`](@ref). """ function embed(M::SpecialEuclideanInGeneralLinear, p, X) G = M.manifold @@ -687,7 +736,11 @@ function embed(M::SpecialEuclideanInGeneralLinear, p, X) Y = allocate_result(G, screw_matrix, nX, hX) nY, hY = submanifold_components(G, Y) copyto!(hY, hX) - copyto!(nY, hp' * nX) + if vector_representation(M.manifold) isa LeftInvariantRepresentation + copyto!(nY, nX) + else + copyto!(nY, hp' * nX) + end @inbounds _padvector!(G, Y) return Y end @@ -699,115 +752,9 @@ 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 ArrayPartition(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 ArrayPartition(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 - -### Special methods for better performance of selected operations - -function exp(M::SpecialEuclidean, p::ArrayPartition, X::ArrayPartition) - M1, M2 = M.manifold.manifolds - return ArrayPartition( - exp(M1.manifold, p.x[1], X.x[1]), - exp(M2.manifold, p.x[2], X.x[2]), - ) -end -function log(M::SpecialEuclidean, p::ArrayPartition, q::ArrayPartition) - M1, M2 = M.manifold.manifolds - return ArrayPartition( - log(M1.manifold, p.x[1], q.x[1]), - log(M2.manifold, p.x[2], q.x[2]), - ) -end -function vee(M::SpecialEuclidean, p::ArrayPartition, X::ArrayPartition) - M1, M2 = M.manifold.manifolds - return vcat(vee(M1.manifold, p.x[1], X.x[1]), vee(M2.manifold, p.x[2], X.x[2])) -end -function get_coordinates( - M::SpecialEuclidean, - p::ArrayPartition, - X::ArrayPartition, - basis::DefaultOrthogonalBasis, -) - M1, M2 = M.manifold.manifolds - return vcat( - get_coordinates(M1.manifold, p.x[1], X.x[1], basis), - get_coordinates(M2.manifold, p.x[2], X.x[2], basis), - ) -end -function hat( - M::SpecialEuclidean{TypeParameter{Tuple{2}}}, - p::ArrayPartition, - c::AbstractVector, -) - M1, M2 = M.manifold.manifolds - return ArrayPartition( - get_vector_orthogonal(M1.manifold, p.x[1], c[SOneTo(2)], ℝ), - get_vector_orthogonal(M2.manifold, p.x[2], c[SA[3]], ℝ), - ) -end -function get_vector( - M::SpecialEuclidean{TypeParameter{Tuple{2}}}, - p::ArrayPartition, - c::AbstractVector, - basis::DefaultOrthogonalBasis, -) - return ArrayPartition( - get_vector(M.manifold.manifolds[1].manifold, p.x[1], c[SOneTo(2)], basis), - get_vector(M.manifold.manifolds[2].manifold, p.x[2], c[SA[3]], basis), - ) -end - -function hat( - M::SpecialEuclidean{TypeParameter{Tuple{3}}}, - p::ArrayPartition, - c::AbstractVector, -) - M1, M2 = M.manifold.manifolds - return ArrayPartition( - get_vector_orthogonal(M1.manifold, p.x[1], c[SOneTo(3)], ℝ), - get_vector_orthogonal(M2.manifold, p.x[2], c[SA[4, 5, 6]], ℝ), - ) -end -function get_vector( - M::SpecialEuclidean{TypeParameter{Tuple{3}}}, - p::ArrayPartition, - c::AbstractVector, - basis::DefaultOrthogonalBasis, -) - return ArrayPartition( - get_vector(M.manifold.manifolds[1].manifold, p.x[1], c[SOneTo(3)], basis), - get_vector(M.manifold.manifolds[2].manifold, p.x[2], c[SA[4, 5, 6]], basis), - ) -end -function compose(::SpecialEuclidean, p::ArrayPartition, q::ArrayPartition) - return ArrayPartition(p.x[2] * q.x[1] + p.x[1], p.x[2] * q.x[2]) -end diff --git a/src/groups/special_linear.jl b/src/groups/special_linear.jl index c29907f080..37462d9a43 100644 --- a/src/groups/special_linear.jl +++ b/src/groups/special_linear.jl @@ -26,7 +26,7 @@ end @inline function active_traits(f, ::SpecialLinear, args...) return merge_traits( - IsGroupManifold(MultiplicationOperation()), + IsGroupManifold(MultiplicationOperation(), LeftInvariantRepresentation()), IsEmbeddedSubmanifold(), HasLeftInvariantMetric(), IsDefaultMetric(EuclideanMetric()), @@ -81,6 +81,9 @@ end inverse_translate_diff(::SpecialLinear, p, q, X, ::LeftForwardAction) = X inverse_translate_diff(::SpecialLinear, p, q, X, ::RightBackwardAction) = p * X / p +adjoint_action!(::SpecialLinear, Y, p, X, ::LeftAction) = copyto!(Y, (p * X) / p) +adjoint_action!(::SpecialLinear, Y, p, X, ::RightAction) = copyto!(Y, p \ X * p) + function inverse_translate_diff!(G::SpecialLinear, Y, p, q, X, conv::ActionDirectionAndSide) return copyto!(Y, inverse_translate_diff(G, p, q, X, conv)) end @@ -152,10 +155,3 @@ function Base.show(io::IO, M::SpecialLinear{Tuple{Int},𝔽}) where {𝔽} n = get_parameter(M.size)[1] return print(io, "SpecialLinear($n, $(𝔽); parameter=:field)") end - -translate_diff(::SpecialLinear, p, q, X, ::LeftForwardAction) = X -translate_diff(::SpecialLinear, p, q, X, ::RightBackwardAction) = p \ X * p - -function translate_diff!(G::SpecialLinear, Y, p, q, X, conv::ActionDirectionAndSide) - return copyto!(Y, translate_diff(G, p, q, X, conv)) -end diff --git a/src/groups/translation_group.jl b/src/groups/translation_group.jl index ff358369e7..5ac303ddbb 100644 --- a/src/groups/translation_group.jl +++ b/src/groups/translation_group.jl @@ -1,5 +1,5 @@ @doc raw""" - TranslationGroup{T,𝔽} <: GroupManifold{Euclidean{T,𝔽},AdditionOperation} + TranslationGroup{T,𝔽} <: GroupManifold{Euclidean{T,𝔽},AdditionOperation,LeftInvariantRepresentation} Translation group ``\mathrm{T}(n)`` represented by translation arrays. @@ -9,23 +9,25 @@ Translation group ``\mathrm{T}(n)`` represented by translation arrays. Generate the translation group on ``𝔽^{n₁,…,nᵢ}`` = `Euclidean(n₁,...,nᵢ; field=𝔽)`, which is isomorphic to the group itself. """ -const TranslationGroup{T,𝔽} = GroupManifold{𝔽,Euclidean{T,𝔽},AdditionOperation} +const TranslationGroup{T,𝔽} = + GroupManifold{𝔽,Euclidean{T,𝔽},AdditionOperation,LeftInvariantRepresentation} function TranslationGroup(n::Int...; field::AbstractNumbers=ℝ, parameter::Symbol=:type) size = wrap_type_parameter(parameter, n) return TranslationGroup{typeof(size),field}( Euclidean(n...; field=field, parameter=parameter), AdditionOperation(), + LeftInvariantRepresentation(), ) end @inline function active_traits(f, M::TranslationGroup, args...) if is_metric_function(f) #pass to Euclidean by default - but keep Group Decorator for the retraction - return merge_traits(IsGroupManifold(M.op), IsExplicitDecorator()) + return merge_traits(IsGroupManifold(M.op, M.vectors), IsExplicitDecorator()) else return merge_traits( - IsGroupManifold(M.op), + IsGroupManifold(M.op, M.vectors), IsDefaultMetric(EuclideanMetric()), active_traits(f, M.manifold, args...), IsExplicitDecorator(), #pass to Euclidean by default/last fallback diff --git a/src/manifolds/FiberBundle.jl b/src/manifolds/FiberBundle.jl index adf08a1d33..fabc1f0c86 100644 --- a/src/manifolds/FiberBundle.jl +++ b/src/manifolds/FiberBundle.jl @@ -9,12 +9,12 @@ Vector transport type on [`FiberBundle`](@ref). # Fields -* `method_horizonal` – vector transport method of the horizontal part (related to manifold M) +* `method_horizontal` – vector transport method of the horizontal part (related to manifold M) * `method_vertical` – vector transport method of the vertical part (related to fibers). The vector transport is derived as a product manifold-style vector transport. The considered product manifold is the product between the manifold ``\mathcal M`` -and the topological vector space isometric to the fiber. +and the space corresponding to the fiber. # Constructor @@ -37,7 +37,7 @@ struct FiberBundleProductVectorTransport{ end function FiberBundleProductVectorTransport( M::AbstractManifold=ManifoldsBase.DefaultManifold(), - fiber::FiberType=ManifoldsBase.TangentSpaceType(); + fiber::FiberType=TangentSpaceType(); vector_transport_method_horizontal::AbstractVectorTransportMethod=default_vector_transport_method( M, ), @@ -370,97 +370,6 @@ function zero_vector!(B::FiberBundle, X, p) return X end -@inline function allocate_result(M::FiberBundle, f::TF) where {TF} - p = allocate_result(M.manifold, f) - X = allocate_result(Fiber(M.manifold, p, M.type), f) - return ArrayPartition(p, X) -end - -function get_vector(M::FiberBundle, p, X, B::AbstractBasis) - n = manifold_dimension(M.manifold) - xp1, xp2 = submanifold_components(M, p) - F = Fiber(M.manifold, xp1, M.type) - return ArrayPartition( - get_vector(M.manifold, xp1, X[1:n], B), - get_vector(F, xp2, X[(n + 1):end], B), - ) -end -function get_vector( - M::FiberBundle, - p, - X, - B::CachedBasis{𝔽,<:AbstractBasis{𝔽},<:FiberBundleBasisData}, -) where {𝔽} - n = manifold_dimension(M.manifold) - xp1, xp2 = submanifold_components(M, p) - F = Fiber(M.manifold, xp1, M.type) - return ArrayPartition( - get_vector(M.manifold, xp1, X[1:n], B.data.base_basis), - get_vector(F, xp2, X[(n + 1):end], B.data.fiber_basis), - ) -end - -function get_vectors( - M::FiberBundle, - p::ArrayPartition, - B::CachedBasis{𝔽,<:AbstractBasis{𝔽},<:FiberBundleBasisData}, -) where {𝔽} - xp1, xp2 = submanifold_components(M, p) - zero_m = zero_vector(M.manifold, xp1) - F = Fiber(M.manifold, xp1, M.type) - zero_f = zero_vector(F, xp1) - vs = typeof(ArrayPartition(zero_m, zero_f))[] - for bv in get_vectors(M.manifold, xp1, B.data.base_basis) - push!(vs, ArrayPartition(bv, zero_f)) - end - for bv in get_vectors(F, xp2, B.data.fiber_basis) - push!(vs, ArrayPartition(zero_m, bv)) - end - return vs -end - -""" - getindex(p::ArrayPartition, M::FiberBundle, s::Symbol) - p[M::FiberBundle, s] - -Access the element(s) at index `s` of a point `p` on a [`FiberBundle`](@ref) `M` by -using the symbols `:point` and `:vector` or `:fiber` for the base and vector or fiber -component, respectively. -""" -@inline function Base.getindex(p::ArrayPartition, M::FiberBundle, s::Symbol) - (s === :point) && return p.x[1] - (s === :vector || s === :fiber) && return p.x[2] - return throw(DomainError(s, "unknown component $s on $M.")) -end - -""" - setindex!(p::ArrayPartition, val, M::FiberBundle, s::Symbol) - p[M::VectorBundle, s] = val - -Set the element(s) at index `s` of a point `p` on a [`FiberBundle`](@ref) `M` to `val` by -using the symbols `:point` and `:fiber` or `:vector` for the base and fiber or vector -component, respectively. - -!!! note - - The *content* of element of `p` is replaced, not the element itself. -""" -@inline function Base.setindex!(x::ArrayPartition, val, M::FiberBundle, s::Symbol) - if s === :point - return copyto!(x.x[1], val) - elseif s === :vector || s === :fiber - return copyto!(x.x[2], val) - else - throw(DomainError(s, "unknown component $s on $M.")) - end -end - function Base.show(io::IO, B::FiberBundle) return print(io, "FiberBundle($(B.type), $(B.manifold), $(B.vector_transport))") end - -@inline function Base.view(x::ArrayPartition, M::FiberBundle, s::Symbol) - (s === :point) && return x.x[1] - (s === :vector || s === :fiber) && return x.x[2] - throw(DomainError(s, "unknown component $s on $M.")) -end diff --git a/src/manifolds/GrassmannStiefel.jl b/src/manifolds/GrassmannStiefel.jl index 14ded44946..934b893b0b 100644 --- a/src/manifolds/GrassmannStiefel.jl +++ b/src/manifolds/GrassmannStiefel.jl @@ -404,27 +404,6 @@ end Base.show(io::IO, p::StiefelPoint) = print(io, "StiefelPoint($(p.value))") Base.show(io::IO, X::StiefelTVector) = print(io, "StiefelTVector($(X.value))") -""" - uniform_distribution(M::Grassmann{<:Any,ℝ}, p) - -Uniform distribution on given (real-valued) [`Grassmann`](@ref) `M`. -Specifically, this is the normalized Haar measure on `M`. -Generated points will be of similar type as `p`. - -The implementation is based on Section 2.5.1 in [Chikuse:2003](@cite); -see also Theorem 2.2.2(iii) in [Chikuse:2003](@cite). -""" -function uniform_distribution(M::Grassmann{<:Any,ℝ}, p) - n, k = get_parameter(M.size) - μ = Distributions.Zeros(n, k) - σ = one(eltype(p)) - Σ1 = Distributions.PDMats.ScalMat(n, σ) - Σ2 = Distributions.PDMats.ScalMat(k, σ) - d = MatrixNormal(μ, Σ1, Σ2) - - return ProjectedPointDistribution(M, d, (M, q, p) -> (q .= svd(p).U), p) -end - # switch order and not dispatch on the _to variant function vector_transport_direction(M::Grassmann, p, X, Y, ::ParallelTransport) return parallel_transport_direction(M, p, X, Y) diff --git a/src/manifolds/MultinomialDoublyStochastic.jl b/src/manifolds/MultinomialDoublyStochastic.jl index 14d09744ca..822c93ddad 100644 --- a/src/manifolds/MultinomialDoublyStochastic.jl +++ b/src/manifolds/MultinomialDoublyStochastic.jl @@ -201,8 +201,8 @@ function project!( ::AbstractMultinomialDoublyStochastic, q, p; - maxiter=100, - tolerance=eps(eltype(p)), + maxiter::Int=100, + tolerance::Real=eps(eltype(p)), ) any(p .<= 0) && throw( DomainError( diff --git a/src/manifolds/MultinomialSymmetricPositiveDefinite.jl b/src/manifolds/MultinomialSymmetricPositiveDefinite.jl index faf8198650..ac4ce954d2 100644 --- a/src/manifolds/MultinomialSymmetricPositiveDefinite.jl +++ b/src/manifolds/MultinomialSymmetricPositiveDefinite.jl @@ -95,19 +95,25 @@ function Random.rand!( p::AbstractMatrix, ) n = get_parameter(M.size)[1] - L = sort(exp.(randn(rng, n))) - V = reduce(hcat, map(xi -> [xi^k for k in 0:(n - 1)], L))' - @static if VERSION < v"1.7" - Vlu = lu(V, Val(false)) - else - Vlu = lu(V, LinearAlgebra.RowNonZero()) + is_spd = false + while !is_spd + L = sort(exp.(randn(rng, n))) + V = reduce(hcat, map(xi -> [xi^k for k in 0:(n - 1)], L))' + @static if VERSION < v"1.7" # COV_EXCL_LINE + Vlu = lu(V, Val(false)) + else + Vlu = lu(V, LinearAlgebra.RowNonZero()) + end + dm = Diagonal(Vlu.U) + uutd = dm \ Vlu.U + random_totally_positive = uutd * dm * Vlu.L + MMDS = MultinomialDoubleStochastic(n) + ds = project(MMDS, random_totally_positive; maxiter=1000) + p .= (ds .+ ds') ./ 2 + if eigmin(p) > 0 + is_spd = true + end end - dm = Diagonal(Vlu.U) - uutd = dm \ Vlu.U - random_totally_positive = uutd * dm * Vlu.L - MMDS = MultinomialDoubleStochastic(n) - ds = project(MMDS, random_totally_positive) - p .= (ds .+ ds') ./ 2 return p end diff --git a/src/manifolds/PowerManifold.jl b/src/manifolds/PowerManifold.jl index 507797386b..7c52066fbb 100644 --- a/src/manifolds/PowerManifold.jl +++ b/src/manifolds/PowerManifold.jl @@ -30,32 +30,6 @@ function PowerManifold( return PowerManifold{𝔽,typeof(M),typeof(size_w),ArrayPowerRepresentation}(M, size_w) end -""" - PowerPointDistribution(M::AbstractPowerManifold, distribution) - -Power distribution on manifold `M`, based on `distribution`. -""" -struct PowerPointDistribution{TM<:AbstractPowerManifold,TD<:MPointDistribution,TX} <: - MPointDistribution{TM} - manifold::TM - distribution::TD - point::TX -end - -""" - PowerFVectorDistribution([type::VectorSpaceFiber], [x], distr) - -Generates a random vector at a `point` from vector space (a fiber of a tangent -bundle) of type `type` using the power distribution of `distr`. - -Vector space type and `point` can be automatically inferred from distribution `distr`. -""" -struct PowerFVectorDistribution{TSpace<:VectorSpaceFiber,TD<:FVectorDistribution} <: - FVectorDistribution{TSpace} - type::TSpace - distribution::TD -end - const PowerManifoldMultidimensional = AbstractPowerManifold{𝔽,<:AbstractManifold{𝔽},ArrayPowerRepresentation} where {𝔽} @@ -64,12 +38,6 @@ Base.:^(M::AbstractManifold, n) = PowerManifold(M, n...) function allocate(::PowerManifoldNestedReplacing, x::AbstractArray{<:SArray}) return similar(x) end -function allocate( - ::PowerManifoldNestedReplacing, - x::AbstractArray{<:ArrayPartition{T,<:NTuple{N,SArray}}}, -) where {T,N} - return similar(x) -end for PowerRepr in [PowerManifoldNested, PowerManifoldNestedReplacing] @eval begin @@ -140,39 +108,6 @@ function manifold_volume(M::PowerManifold) return manifold_volume(M.manifold)^prod(size) end -function Random.rand(rng::AbstractRNG, d::PowerFVectorDistribution) - fv = zero_vector(d.type.manifold, d.type.point) - Distributions._rand!(rng, d, fv) - return fv -end -function Random.rand(rng::AbstractRNG, d::PowerPointDistribution) - x = allocate_result(d.manifold, rand, d.point) - Distributions._rand!(rng, d, x) - return x -end - -function Distributions._rand!( - rng::AbstractRNG, - d::PowerFVectorDistribution, - v::AbstractArray, -) - PM = d.type.manifold - rep_size = representation_size(PM.manifold) - for i in get_iterator(d.type.manifold) - copyto!(d.distribution.type.point, _read(PM, rep_size, d.type.point, i)) - Distributions._rand!(rng, d.distribution, _read(PM, rep_size, v, i)) - end - return v -end -function Distributions._rand!(rng::AbstractRNG, d::PowerPointDistribution, x::AbstractArray) - M = d.manifold - rep_size = representation_size(M.manifold) - for i in get_iterator(M) - Distributions._rand!(rng, d.distribution, _write(M, rep_size, x, i)) - end - return x -end - Base.@propagate_inbounds @inline function _read( ::PowerManifoldMultidimensional, rep_size::Tuple, @@ -189,23 +124,6 @@ Base.@propagate_inbounds @inline function _read( ) where {N} return x[i...] end -Base.@propagate_inbounds @inline function _read( - ::PowerManifoldMultidimensional, - rep_size::Tuple, - x::HybridArray, - i::Tuple, -) - return x[rep_size_to_colons(rep_size)..., i...] -end -Base.@propagate_inbounds @inline function _read( - ::PowerManifoldMultidimensional, - rep_size::Tuple{}, - x::HybridArray, - i::NTuple{N,Int}, -) where {N} - # disambiguation - return x[i...] -end function Base.view( p::AbstractArray, @@ -287,9 +205,6 @@ function Base.show( return print(io, "PowerManifold($(M.manifold), $(join(size, ", ")))") end -Distributions.support(tvd::PowerFVectorDistribution) = FVectorSupport(tvd.type) -Distributions.support(d::PowerPointDistribution) = MPointSupport(d.manifold) - @doc raw""" volume_density(M::PowerManifold, p, X) diff --git a/src/manifolds/ProductManifold.jl b/src/manifolds/ProductManifold.jl index cc29c2eea5..60ef271bed 100644 --- a/src/manifolds/ProductManifold.jl +++ b/src/manifolds/ProductManifold.jl @@ -7,55 +7,6 @@ function allocate_coordinates(::ProductManifold, p, T, n::Int) return allocate(submanifold_component(p, 1), T, n) end -""" - ProductFVectorDistribution([type::VectorSpaceFiber], [x], distrs...) - -Generates a random vector at point `x` from vector space (a fiber of a tangent -bundle) of type `type` using the product distribution of given distributions. - -Vector space type and `x` can be automatically inferred from distributions `distrs`. -""" -struct ProductFVectorDistribution{ - TSpace<:VectorSpaceFiber{<:Any,<:ProductManifold}, - TD<:(NTuple{N,Distribution} where {N}), -} <: FVectorDistribution{TSpace} - type::TSpace - distributions::TD -end - -""" - ProductPointDistribution(M::ProductManifold, distributions) - -Product distribution on manifold `M`, combined from `distributions`. -""" -struct ProductPointDistribution{ - TM<:ProductManifold, - TD<:(NTuple{N,Distribution} where {N}), -} <: MPointDistribution{TM} - manifold::TM - distributions::TD -end - -function adjoint_Jacobi_field( - M::ProductManifold, - p::ArrayPartition, - q::ArrayPartition, - t, - X::ArrayPartition, - β::Tβ, -) where {Tβ} - return ArrayPartition( - map( - adjoint_Jacobi_field, - M.manifolds, - submanifold_components(M, p), - submanifold_components(M, q), - ntuple(_ -> t, length(M.manifolds)), - submanifold_components(M, X), - ntuple(_ -> β, length(M.manifolds)), - )..., - ) -end function adjoint_Jacobi_field!(M::ProductManifold, Y, p, q, t, X, β::Tβ) where {Tβ} map( adjoint_Jacobi_field!, @@ -80,26 +31,6 @@ to the corresponding entry in `p`) separately. """ flat(::ProductManifold, ::Any...) -function jacobi_field( - M::ProductManifold, - p::ArrayPartition, - q::ArrayPartition, - t, - X::ArrayPartition, - β::Tβ, -) where {Tβ} - return ArrayPartition( - map( - jacobi_field, - M.manifolds, - submanifold_components(M, p), - submanifold_components(M, q), - ntuple(_ -> t, length(M.manifolds)), - submanifold_components(M, X), - ntuple(_ -> β, length(M.manifolds)), - )..., - ) -end function jacobi_field!(M::ProductManifold, Y, p, q, t, X, β::Tβ) where {Tβ} map( jacobi_field!, @@ -122,73 +53,6 @@ Return the volume of [`ProductManifold`](@extref `ManifoldsBase.ProductManifold` """ manifold_volume(M::ProductManifold) = mapreduce(manifold_volume, *, M.manifolds) -function ProductFVectorDistribution(distributions::FVectorDistribution...) - M = ProductManifold(map(d -> support(d).space.manifold, distributions)...) - fiber_type = support(distributions[1]).space.fiber_type - if !all(d -> support(d).space.fiber_type == fiber_type, distributions) - error( - "Not all distributions have support in vector spaces of the same type, which is currently not supported", - ) - end - # Probably worth considering sum spaces in the future? - p = ArrayPartition(map(d -> support(d).space.point, distributions)...) - return ProductFVectorDistribution(Fiber(M, p, fiber_type), distributions) -end - -function ProductPointDistribution(M::ProductManifold, distributions::MPointDistribution...) - return ProductPointDistribution{typeof(M),typeof(distributions)}(M, distributions) -end -function ProductPointDistribution(distributions::MPointDistribution...) - M = ProductManifold(map(d -> support(d).manifold, distributions)...) - return ProductPointDistribution(M, distributions...) -end - -function Random.rand(rng::AbstractRNG, d::ProductPointDistribution) - return ArrayPartition(map(d -> rand(rng, d), d.distributions)...) -end -function Random.rand(rng::AbstractRNG, d::ProductFVectorDistribution) - return ArrayPartition(map(d -> rand(rng, d), d.distributions)...) -end - -function Distributions._rand!( - rng::AbstractRNG, - d::ProductPointDistribution, - x::AbstractArray{<:Number}, -) - return copyto!(x, rand(rng, d)) -end -function Distributions._rand!( - rng::AbstractRNG, - d::ProductPointDistribution, - p::ArrayPartition, -) - map( - (t1, t2) -> Distributions._rand!(rng, t1, t2), - d.distributions, - submanifold_components(d.manifold, p), - ) - return p -end -function Distributions._rand!( - rng::AbstractRNG, - d::ProductFVectorDistribution, - v::AbstractArray{<:Number}, -) - return copyto!(v, rand(rng, d)) -end -function Distributions._rand!( - rng::AbstractRNG, - d::ProductFVectorDistribution, - X::ArrayPartition, -) - map( - t -> Distributions._rand!(rng, t[1], t[2]), - d.distributions, - submanifold_components(d.space.manifold, X), - ) - return X -end - @doc raw""" Y = riemannian_Hessian(M::ProductManifold, p, G, H, X) riemannian_Hessian!(M::ProductManifold, Y, p, G, H, X) @@ -224,21 +88,6 @@ separately """ sharp(::ProductManifold, ::Any...) -Distributions.support(d::ProductPointDistribution) = MPointSupport(d.manifold) -function Distributions.support(tvd::ProductFVectorDistribution) - return FVectorSupport(tvd.type) -end - -function uniform_distribution(M::ProductManifold) - return ProductPointDistribution(M, map(uniform_distribution, M.manifolds)) -end -function uniform_distribution(M::ProductManifold, p) - return ProductPointDistribution( - M, - map(uniform_distribution, M.manifolds, submanifold_components(M, p)), - ) -end - @doc raw""" volume_density(M::ProductManifold, p, X) diff --git a/src/manifolds/ProjectiveSpace.jl b/src/manifolds/ProjectiveSpace.jl index 2855deb948..6d3ea67d4a 100644 --- a/src/manifolds/ProjectiveSpace.jl +++ b/src/manifolds/ProjectiveSpace.jl @@ -506,17 +506,6 @@ function Base.show(io::IO, M::ArrayProjectiveSpace{<:Tuple,𝔽}) where {𝔽} return print(io, "ArrayProjectiveSpace($(join(n, ", ")); field=$(𝔽), parameter=:field)") end -""" - uniform_distribution(M::ProjectiveSpace{<:Any,ℝ}, p) - -Uniform distribution on given [`ProjectiveSpace`](@ref) `M`. Generated points will be of -similar type as `p`. -""" -function uniform_distribution(M::ProjectiveSpace{<:Any,ℝ}, p) - d = Distributions.MvNormal(zero(p), 1.0 * I) - return ProjectedPointDistribution(M, d, project!, p) -end - @doc raw""" parallel_transport_to(M::AbstractProjectiveSpace, p, X, q) diff --git a/src/manifolds/Rotations.jl b/src/manifolds/Rotations.jl index fccece31c5..204a23faf6 100644 --- a/src/manifolds/Rotations.jl +++ b/src/manifolds/Rotations.jl @@ -17,29 +17,6 @@ function Rotations(n::Int; parameter::Symbol=:type) return Rotations{typeof(size)}(size) end -""" - NormalRotationDistribution(M::Rotations, d::Distribution, x::TResult) - -Distribution that returns a random point on the manifold [`Rotations`](@ref) -`M`. Random point is generated using base distribution `d` and the type -of the result is adjusted to `TResult`. - -See [`normal_rotation_distribution`](@ref) for details. -""" -struct NormalRotationDistribution{TResult,TM<:Rotations,TD<:Distribution} <: - MPointDistribution{TM} - manifold::TM - distr::TD -end - -function NormalRotationDistribution( - M::Rotations, - d::Distribution, - x::TResult, -) where {TResult} - return NormalRotationDistribution{TResult,typeof(M),typeof(d)}(M, d) -end - @doc raw""" angles_4d_skew_sym_matrix(A) @@ -259,11 +236,7 @@ and second columns are swapped. The argument `p` is used to determine the type of returned points. """ -function normal_rotation_distribution(M::Rotations, p, σ::Real) - n = get_parameter(M.size)[1] - d = Distributions.MvNormal(zeros(n * n), σ * I) - return NormalRotationDistribution(M, d, p) -end +function normal_rotation_distribution end @doc raw""" project(M::Rotations, p; check_det = true) @@ -297,19 +270,6 @@ function project!(M::Rotations, q, p; check_det::Bool=true) return q end -function Random.rand( - rng::AbstractRNG, - d::NormalRotationDistribution{TResult,<:Rotations}, -) where {TResult} - n = get_parameter(d.manifold.size)[1] - return if n == 1 - convert(TResult, ones(1, 1)) - else - A = reshape(rand(rng, d.distr), (n, n)) - convert(TResult, _fix_random_rotation(A)) - end -end - function Random.rand!( rng::AbstractRNG, M::Rotations, @@ -336,24 +296,6 @@ function Random.rand!( return pX end -function Distributions._rand!( - rng::AbstractRNG, - d::NormalRotationDistribution, - x::AbstractArray{<:Real}, -) - return copyto!(x, rand(rng, d)) -end - -function _fix_random_rotation(A::AbstractMatrix) - s = diag(sign.(qr(A).R)) - D = Diagonal(s) - C = qr(A).Q * D - if det(C) < 0 - C[:, [1, 2]] = C[:, [2, 1]] - end - return C -end - @doc raw""" parallel_transport_direction(M::Rotations, p, X, d) @@ -454,8 +396,6 @@ function sectional_curvature_min(::Rotations) return 0.0 end -Distributions.support(d::NormalRotationDistribution) = MPointSupport(d.manifold) - @doc raw""" Weingarten(M::Rotations, p, X, V) diff --git a/src/manifolds/Sphere.jl b/src/manifolds/Sphere.jl index a9981cd1b3..9cc8948d0c 100644 --- a/src/manifolds/Sphere.jl +++ b/src/manifolds/Sphere.jl @@ -518,18 +518,6 @@ function Base.show(io::IO, M::ArraySphere{<:Tuple,𝔽}) where {𝔽} return print(io, "ArraySphere($(join(n, ", ")); field=$(𝔽), parameter=:field)") end -""" - uniform_distribution(M::Sphere{n,ℝ}, p) where {n} - -Uniform distribution on given [`Sphere`](@ref) `M`. Generated points will be of -similar type as `p`. -""" -function uniform_distribution(M::Sphere{<:Any,ℝ}, p) - n = get_parameter(M.size)[1] - d = Distributions.MvNormal(zero(p), 1.0 * I) - return ProjectedPointDistribution(M, d, project!, p) -end - @doc raw""" parallel_transport_to(M::AbstractSphere, p, X, q) diff --git a/src/manifolds/Stiefel.jl b/src/manifolds/Stiefel.jl index b8243dcc26..0f72f2ce0e 100644 --- a/src/manifolds/Stiefel.jl +++ b/src/manifolds/Stiefel.jl @@ -522,27 +522,6 @@ function Base.show(io::IO, M::Stiefel{Tuple{Int,Int},𝔽}) where {𝔽} return print(io, "Stiefel($(n), $(k), $(𝔽); parameter=:field)") end -""" - uniform_distribution(M::Stiefel{<:Any,ℝ}, p) - -Uniform distribution on given (real-valued) [`Stiefel`](@ref) `M`. -Specifically, this is the normalized Haar and Hausdorff measure on `M`. -Generated points will be of similar type as `p`. - -The implementation is based on Section 2.5.1 in [Chikuse:2003](@cite); -see also Theorem 2.2.1(iii) in [Chikuse:2003](@cite). -""" -function uniform_distribution(M::Stiefel{<:Any,ℝ}, p) - n, k = get_parameter(M.size) - μ = Distributions.Zeros(n, k) - σ = one(eltype(p)) - Σ1 = Distributions.PDMats.ScalMat(n, σ) - Σ2 = Distributions.PDMats.ScalMat(k, σ) - d = MatrixNormal(μ, Σ1, Σ2) - - return ProjectedPointDistribution(M, d, project!, p) -end - @doc raw""" vector_transport_direction(::Stiefel, p, X, d, ::DifferentiatedRetractionVectorTransport{CayleyRetraction}) diff --git a/src/manifolds/Symplectic.jl b/src/manifolds/Symplectic.jl index 944ee44218..65ffce5312 100644 --- a/src/manifolds/Symplectic.jl +++ b/src/manifolds/Symplectic.jl @@ -696,7 +696,7 @@ function project_normal!( end @doc raw""" - rand(::SymplecticStiefel; vector_at=nothing, σ=1.0) + rand(::SymplecticStiefel; vector_at=nothing, σ::Real=1.0) Generate a random point on ``\mathrm{Sp}(2n)`` or a random tangent vector ``X \in T_p\mathrm{Sp}(2n)`` if `vector_at` is set to @@ -718,20 +718,15 @@ and then symmetrizes it as `S = S + S'`. Then ``S`` is normalized to have Frobenius norm of `σ` and `X = pJS` is returned, where `J` is the [`SymplecticElement`](@ref). """ -rand(SymplecticMatrices; σ::Rieal=1.0, kwargs...) +rand(SymplecticMatrices; σ::Real=1.0, kwargs...) function Random.rand!( rng::AbstractRNG, M::SymplecticMatrices, pX; vector_at=nothing, - hamiltonian_norm=nothing, - σ=hamiltonian_norm === nothing ? 1.0 : hamiltonian_norm, + σ::Real=1.0, ) - !(hamiltonian_norm === nothing) && Base.depwarn( - Random.rand!, - "hamiltonian_norm is deprecated as a keyword, please use the default σ.", - ) n = get_parameter(M.size)[1] if vector_at === nothing rand!(rng, HamiltonianMatrices(2n), pX; σ=σ) @@ -743,7 +738,7 @@ function Random.rand!( end end -function random_vector!(M::SymplecticMatrices, X, p; σ=1.0) +function random_vector!(M::SymplecticMatrices, X, p; σ::Real=1.0) n = get_parameter(M.size)[1] # Generate random symmetric matrix: randn!(X) diff --git a/src/manifolds/SymplecticGrassmann.jl b/src/manifolds/SymplecticGrassmann.jl index c5051af2e2..92b6da4523 100644 --- a/src/manifolds/SymplecticGrassmann.jl +++ b/src/manifolds/SymplecticGrassmann.jl @@ -72,14 +72,9 @@ struct SymplecticGrassmann{T,𝔽} <: AbstractDecoratorManifold{𝔽} size::T end -function SymplecticGrassmann( - two_n::Int, - two_k::Int, - field::AbstractNumbers=ℝ; - parameter::Symbol=:type, -) +function SymplecticGrassmann(two_n::Int, two_k::Int; parameter::Symbol=:type) size = wrap_type_parameter(parameter, (div(two_n, 2), div(two_k, 2))) - return SymplecticGrassmann{typeof(size),field}(size) + return SymplecticGrassmann{typeof(size),ℝ}(size) end function active_traits(f, ::SymplecticGrassmann, args...) @@ -105,12 +100,12 @@ function manifold_dimension(M::SymplecticGrassmann{<:Any,ℝ}) return 4 * (n - k) * k end -function Base.show(io::IO, ::SymplecticGrassmann{TypeParameter{Tuple{n,k}},𝔽}) where {n,k,𝔽} - return print(io, "SymplecticGrassmann($(2n), $(2k); field=$(𝔽))") +function Base.show(io::IO, ::SymplecticGrassmann{TypeParameter{Tuple{n,k}}}) where {n,k} + return print(io, "SymplecticGrassmann($(2n), $(2k))") end -function Base.show(io::IO, M::SymplecticGrassmann{Tuple{Int,Int},𝔽}) where {𝔽} +function Base.show(io::IO, M::SymplecticGrassmann{Tuple{Int,Int}}) n, k = get_parameter(M.size) - return print(io, "SymplecticGrassmann($(2n), $(2k); field=$(𝔽); parameter=:field)") + return print(io, "SymplecticGrassmann($(2n), $(2k); parameter=:field)") end # diff --git a/src/manifolds/SymplecticGrassmannProjector.jl b/src/manifolds/SymplecticGrassmannProjector.jl index 4a623f783a..22726a010b 100644 --- a/src/manifolds/SymplecticGrassmannProjector.jl +++ b/src/manifolds/SymplecticGrassmannProjector.jl @@ -80,12 +80,12 @@ embed(::SymplecticGrassmann, p::ProjectorPoint) = p.value embed(::SymplecticGrassmann, p::ProjectorPoint, X::ProjectorTVector) = X.value function get_embedding( - ::SymplecticGrassmann{TypeParameter{Tuple{n,k}},𝔽}, + ::SymplecticGrassmann{TypeParameter{Tuple{n,k}}}, p::ProjectorPoint, -) where {n,k,𝔽} - return Euclidean(2n, 2n; field=𝔽) +) where {n,k} + return Euclidean(2n, 2n;) end -function get_embedding(M::SymplecticGrassmann{Tuple{Int,Int},𝔽}, ::ProjectorPoint) where {𝔽} +function get_embedding(M::SymplecticGrassmann{Tuple{Int,Int}}, ::ProjectorPoint) n, _ = get_parameter(M.size) - return Euclidean(2n, 2n; field=𝔽, parameter=:field) + return Euclidean(2n, 2n; parameter=:field) end diff --git a/src/manifolds/SymplecticGrassmannStiefel.jl b/src/manifolds/SymplecticGrassmannStiefel.jl index b1a52d25f2..0bf394802d 100644 --- a/src/manifolds/SymplecticGrassmannStiefel.jl +++ b/src/manifolds/SymplecticGrassmannStiefel.jl @@ -54,12 +54,12 @@ function exp!(M::SymplecticGrassmann, q, p, X) return q end -function get_embedding(::SymplecticGrassmann{TypeParameter{Tuple{n,k}},𝔽}) where {n,k,𝔽} - return SymplecticStiefel(2n, 2k, 𝔽) +function get_embedding(::SymplecticGrassmann{TypeParameter{Tuple{n,k}}}) where {n,k} + return SymplecticStiefel(2n, 2k) end -function get_embedding(M::SymplecticGrassmann{Tuple{Int,Int},𝔽}) where {𝔽} +function get_embedding(M::SymplecticGrassmann{Tuple{Int,Int}}) n, k = get_parameter(M.size) - return SymplecticStiefel(2n, 2k, 𝔽; parameter=:field) + return SymplecticStiefel(2n, 2k; parameter=:field) end @doc raw""" diff --git a/src/manifolds/SymplecticStiefel.jl b/src/manifolds/SymplecticStiefel.jl index 46809b01f2..76c50c3437 100644 --- a/src/manifolds/SymplecticStiefel.jl +++ b/src/manifolds/SymplecticStiefel.jl @@ -47,12 +47,7 @@ struct SymplecticStiefel{T,𝔽} <: AbstractDecoratorManifold{𝔽} size::T end -function SymplecticStiefel( - two_n::Int, - two_k::Int, - field::AbstractNumbers=ℝ; - parameter::Symbol=:type, -) +function SymplecticStiefel(two_n::Int, two_k::Int; parameter::Symbol=:type) two_n % 2 == 0 || throw( ArgumentError( "The first matrix size of the symplectic Stiefel manifold must be even, but was $(two_n).", @@ -64,7 +59,7 @@ function SymplecticStiefel( ), ) size = wrap_type_parameter(parameter, (div(two_n, 2), div(two_k, 2))) - return SymplecticStiefel{typeof(size),field}(size) + return SymplecticStiefel{typeof(size),ℝ}(size) end function active_traits(f, ::SymplecticStiefel, args...) @@ -72,7 +67,7 @@ function active_traits(f, ::SymplecticStiefel, args...) end # Define Stiefel as the array fallback -ManifoldsBase.@default_manifold_fallbacks SymplecticStiefel StiefelPoint StiefelTVector value value +ManifoldsBase.@default_manifold_fallbacks SymplecticStiefel{<:Any,ℝ} StiefelPoint StiefelTVector value value function ManifoldsBase.default_inverse_retraction_method(::SymplecticStiefel) return CayleyInverseRetraction() @@ -506,13 +501,8 @@ function Random.rand!( M::SymplecticStiefel, pX; vector_at=nothing, - hamiltonian_norm=nothing, - σ=hamiltonian_norm === nothing ? 1.0 : hamiltonian_norm, + σ::Real=1.0, ) - !(hamiltonian_norm === nothing) && Base.depwarn( - Random.rand!, - "hamiltonian_norm is deprecated as a keyword, please use the default σ.", - ) n, k = get_parameter(M.size) if vector_at === nothing canonical_project!(M, pX, rand(rng, SymplecticMatrices(2n); σ=σ)) @@ -620,12 +610,12 @@ function riemannian_gradient!( return X end -function Base.show(io::IO, ::SymplecticStiefel{TypeParameter{Tuple{n,k}},𝔽}) where {n,k,𝔽} - return print(io, "SymplecticStiefel($(2n), $(2k); field=$(𝔽))") +function Base.show(io::IO, ::SymplecticStiefel{TypeParameter{Tuple{n,k}}}) where {n,k} + return print(io, "SymplecticStiefel($(2n), $(2k))") end -function Base.show(io::IO, M::SymplecticStiefel{Tuple{Int,Int},𝔽}) where {𝔽} +function Base.show(io::IO, M::SymplecticStiefel{Tuple{Int,Int}}) n, k = get_parameter(M.size) - return print(io, "SymplecticStiefel($(2n), $(2k); field=$(𝔽); parameter=:field)") + return print(io, "SymplecticStiefel($(2n), $(2k); parameter=:field)") end @doc raw""" diff --git a/src/manifolds/VectorBundle.jl b/src/manifolds/VectorBundle.jl index 00c5e7995d..70fab52fc0 100644 --- a/src/manifolds/VectorBundle.jl +++ b/src/manifolds/VectorBundle.jl @@ -1,11 +1,4 @@ -""" - const VectorBundleVectorTransport = FiberBundleProductVectorTransport - -Deprecated: an alias for `FiberBundleProductVectorTransport`. -""" -const VectorBundleVectorTransport = FiberBundleProductVectorTransport - """ fiber_bundle_transport(M::AbstractManifold, fiber::FiberType) @@ -383,36 +376,4 @@ function vector_transport_to!( return Y end -function _vector_transport_direction( - M::VectorBundle, - p, - X, - d, - m::FiberBundleProductVectorTransport, -) - px, pVx = submanifold_components(M.manifold, p) - VXM, VXF = submanifold_components(M.manifold, X) - dx, dVx = submanifold_components(M.manifold, d) - return ArrayPartition( - vector_transport_direction(M.manifold, px, VXM, dx, m.method_horizontal), - bundle_transport_tangent_direction(M, px, pVx, VXF, dx, m.method_vertical), - ) -end - -function _vector_transport_to( - M::VectorBundle, - p, - X, - q, - m::FiberBundleProductVectorTransport, -) - px, pVx = submanifold_components(M.manifold, p) - VXM, VXF = submanifold_components(M.manifold, X) - qx, qVx = submanifold_components(M.manifold, q) - return ArrayPartition( - vector_transport_to(M.manifold, px, VXM, qx, m.method_horizontal), - bundle_transport_tangent_to(M, px, pVx, VXF, qx, m.method_vertical), - ) -end - Base.show(io::IO, vb::CotangentBundle) = print(io, "CotangentBundle($(vb.manifold))") diff --git a/src/product_representations.jl b/src/product_representations.jl deleted file mode 100644 index 8b13789179..0000000000 --- a/src/product_representations.jl +++ /dev/null @@ -1 +0,0 @@ - diff --git a/src/statistics.jl b/src/statistics.jl index d2800764ca..3e08cf3780 100644 --- a/src/statistics.jl +++ b/src/statistics.jl @@ -1,9 +1,3 @@ -""" - AbstractEstimationMethod - -Deprecated alias for `AbstractApproximationMethod` -""" -const AbstractEstimationMethod = AbstractApproximationMethod _unit_weights(n::Int) = StatsBase.UnitWeights{Float64}(n) @@ -407,16 +401,8 @@ function Statistics.mean!( x::AbstractVector, w::AbstractVector, e::ExtrinsicEstimation; - extrinsic_method::Union{AbstractEstimationMethod,Nothing}=nothing, kwargs..., ) - if !isnothing(extrinsic_method) - Base.depwarn( - "The Keyword Argument `extrinsic_method` is deprecated use `ExtrinsicEstimators` field instead", - :mean!, - ) - e = ExtrinsicEstimation(extrinsic_method) - end embedded_x = map(p -> embed(M, p), x) embedded_y = mean(get_embedding(M), embedded_x, w, e.extrinsic_estimation; kwargs...) project!(M, y, embedded_y) @@ -685,16 +671,8 @@ function Statistics.median!( x::AbstractVector, w::AbstractVector, e::ExtrinsicEstimation; - extrinsic_method=nothing, kwargs..., ) - if !isnothing(extrinsic_method) - Base.depwarn( - "The Keyword Argument `extrinsic_method` is deprecated use `ExtrinsicEstimators` field instead", - :median!, - ) - e = ExtrinsicEstimation(extrinsic_method) - end embedded_x = map(p -> embed(M, p), x) embedded_y = median(get_embedding(M), embedded_x, w, e.extrinsic_estimation; kwargs...) project!(M, y, embedded_y) diff --git a/test/ambiguities.jl b/test/ambiguities.jl index 90f707714a..fedf91fc10 100644 --- a/test/ambiguities.jl +++ b/test/ambiguities.jl @@ -17,7 +17,7 @@ end # Interims solution until we follow what was proposed in # https://discourse.julialang.org/t/avoid-ambiguities-with-individual-number-element-identity/62465/2 fmbs = filter(x -> !any(has_type_in_signature.(x, Identity)), mbs) - FMBS_LIMIT = 38 + FMBS_LIMIT = 33 println("Number of ManifoldsBase.jl ambiguities: $(length(fmbs))") @test length(fmbs) <= FMBS_LIMIT if length(fmbs) > FMBS_LIMIT @@ -30,7 +30,7 @@ end # Interims solution until we follow what was proposed in # https://discourse.julialang.org/t/avoid-ambiguities-with-individual-number-element-identity/62465/2 fms = filter(x -> !any(has_type_in_signature.(x, Identity)), ms) - FMS_LIMIT = 53 + FMS_LIMIT = 46 println("Number of Manifolds.jl ambiguities: $(length(fms))") if length(fms) > FMS_LIMIT for amb in fms diff --git a/test/groups/general_linear.jl b/test/groups/general_linear.jl index 1af110b040..0f1ffa73a7 100644 --- a/test/groups/general_linear.jl +++ b/test/groups/general_linear.jl @@ -2,6 +2,8 @@ include("../header.jl") include("group_utils.jl") using NLsolve +using Manifolds: LeftForwardAction, RightBackwardAction + @testset "General Linear group" begin @testset "basic properties" begin G = GeneralLinear(3) diff --git a/test/groups/group_operation_action.jl b/test/groups/group_operation_action.jl index ec7e9d569e..fbede58dfe 100644 --- a/test/groups/group_operation_action.jl +++ b/test/groups/group_operation_action.jl @@ -6,7 +6,11 @@ using Manifolds: LeftForwardAction, LeftBackwardAction, RightForwardAction, RightBackwardAction @testset "Group operation action" begin - G = GroupManifold(NotImplementedManifold(), Manifolds.MultiplicationOperation()) + G = GroupManifold( + NotImplementedManifold(), + Manifolds.MultiplicationOperation(), + Manifolds.LeftInvariantRepresentation(), + ) A_left_fwd = GroupOperationAction(G) A_left_back = GroupOperationAction(G, LeftBackwardAction()) A_right_fwd = GroupOperationAction(G, RightForwardAction()) diff --git a/test/groups/group_utils.jl b/test/groups/group_utils.jl index 4254997dfc..6c3c104c6a 100644 --- a/test/groups/group_utils.jl +++ b/test/groups/group_utils.jl @@ -24,7 +24,7 @@ if !s end function active_traits(f, M::DefaultTransparencyGroup, args...) return merge_traits( - Manifolds.IsGroupManifold(M.op), + Manifolds.IsGroupManifold(M.op, LeftInvariantRepresentation()), active_traits(f, M.manifold, args...), ) end diff --git a/test/groups/groups_general.jl b/test/groups/groups_general.jl index 710167904b..64da8572a9 100644 --- a/test/groups/groups_general.jl +++ b/test/groups/groups_general.jl @@ -9,7 +9,11 @@ using Manifolds: @testset "General group tests" begin @testset "Not implemented operation" begin - G = GroupManifold(NotImplementedManifold(), NotImplementedOperation()) + G = GroupManifold( + NotImplementedManifold(), + NotImplementedOperation(), + Manifolds.LeftInvariantRepresentation(), + ) @test repr(G) == "GroupManifold(NotImplementedManifold(), NotImplementedOperation())" p = [1.0, 2.0] @@ -57,8 +61,14 @@ using Manifolds: @test !isapprox(G, eg, Identity(AdditionOperation())) @test !isapprox(G, Identity(AdditionOperation()), eg) - @test NotImplementedOperation(NotImplementedManifold()) === G - @test (NotImplementedOperation())(NotImplementedManifold()) === G + @test NotImplementedOperation( + NotImplementedManifold(), + Manifolds.LeftInvariantRepresentation(), + ) === G + @test (NotImplementedOperation())( + NotImplementedManifold(), + Manifolds.LeftInvariantRepresentation(), + ) === G @test_throws ErrorException hat(Rotations(3), eg, [1, 2, 3]) @test_throws ErrorException hat!(Rotations(3), randn(3, 3), eg, [1, 2, 3]) @@ -80,6 +90,8 @@ using Manifolds: @test_throws MethodError inv!(G, p, p) @test_throws MethodError inv!(G, p, eg) @test_throws MethodError inv(G, p) + @test_throws MethodError inv_diff(G, p, X) + @test_throws MethodError inv_diff!(G, X, p, X) # no function defined to return the identity array representation @test_throws MethodError copyto!(G, p, eg) @@ -108,26 +120,12 @@ using Manifolds: @test_throws MethodError inverse_translate!(G, p, p, p, LeftForwardAction()) @test_throws MethodError inverse_translate!(G, p, p, p, RightBackwardAction()) - @test_throws MethodError translate_diff(G, p, p, X) - @test_throws MethodError translate_diff(G, p, p, X, LeftForwardAction()) - @test_throws MethodError translate_diff(G, p, p, X, RightBackwardAction()) - @test_throws MethodError translate_diff!(G, X, p, p, X) - @test_throws MethodError translate_diff!(G, X, p, p, X, LeftForwardAction()) - @test_throws MethodError translate_diff!(G, X, p, p, X, RightBackwardAction()) - - @test_throws MethodError inverse_translate_diff(G, p, p, X) - @test_throws MethodError inverse_translate_diff(G, p, p, X, LeftForwardAction()) - @test_throws MethodError inverse_translate_diff(G, p, p, X, RightBackwardAction()) - @test_throws MethodError inverse_translate_diff!(G, X, p, p, X) - @test_throws MethodError inverse_translate_diff!(G, X, p, p, X, LeftForwardAction()) - @test_throws MethodError inverse_translate_diff!( - G, - X, - p, - p, - X, - RightBackwardAction(), - ) + @test_throws MethodError adjoint_action(G, p, X) + @test_throws MethodError adjoint_action(G, p, X, LeftAction()) + @test_throws MethodError adjoint_action(G, p, X, RightAction()) + @test_throws MethodError adjoint_action!(G, p, X, X) + @test_throws MethodError adjoint_action!(G, p, X, X, LeftAction()) + @test_throws MethodError adjoint_action!(G, p, X, X, RightAction()) @test_throws MethodError exp_lie(G, X) @test_throws MethodError exp_lie!(G, p, X) @@ -140,7 +138,11 @@ using Manifolds: @test switch_direction(LeftAction()) === RightAction() @test switch_direction(RightAction()) === LeftAction() - G = GroupManifold(NotImplementedManifold(), NotImplementedOperation()) + G = GroupManifold( + NotImplementedManifold(), + NotImplementedOperation(), + Manifolds.LeftInvariantRepresentation(), + ) @test Manifolds._action_order(G, 1, 2, LeftForwardAction()) === (1, 2) @test Manifolds._action_order(G, 1, 2, RightBackwardAction()) === (2, 1) end @@ -151,7 +153,11 @@ using Manifolds: end @testset "Addition operation" begin - G = GroupManifold(NotImplementedManifold(), Manifolds.AdditionOperation()) + G = GroupManifold( + NotImplementedManifold(), + Manifolds.AdditionOperation(), + Manifolds.LeftInvariantRepresentation(), + ) test_group( G, [[1.0, 2.0], [3.0, 4.0], [5.0, 6.0]], @@ -197,7 +203,11 @@ using Manifolds: end @testset "Multiplication operation" begin - G = GroupManifold(NotImplementedManifold(), Manifolds.MultiplicationOperation()) + G = GroupManifold( + NotImplementedManifold(), + Manifolds.MultiplicationOperation(), + Manifolds.LeftInvariantRepresentation(), + ) test_group( G, [[2.0 1.0; 3.0 4.0], [3.0 2.0; 4.0 5.0], [4.0 3.0; 5.0 6.0]], diff --git a/test/groups/product_group.jl b/test/groups/product_group.jl index b1e1f7e978..2fd56f6b41 100644 --- a/test/groups/product_group.jl +++ b/test/groups/product_group.jl @@ -7,8 +7,11 @@ using RecursiveArrayTools Tn = TranslationGroup(2) Rn = Rotations(3) M = ProductManifold(SOn, Tn) - G = ProductGroup(M) - @test_throws ErrorException ProductGroup(ProductManifold(Rotations(3), Stiefel(3, 2))) + G = ProductGroup(M, Manifolds.LeftInvariantRepresentation()) + @test_throws ErrorException ProductGroup( + ProductManifold(Rotations(3), Stiefel(3, 2)), + Manifolds.LeftInvariantRepresentation(), + ) @test G isa ProductGroup @test submanifold(G, 1) === SOn @test submanifold(G, 2) === Tn @@ -30,7 +33,7 @@ using RecursiveArrayTools end pts = [ArrayPartition(tp...) for tp in tuple_pts] - X_pts = [ArrayPartition(tuple_v...)] + X_pts = [ArrayPartition(tuple_v...), ArrayPartition(tuple_v...)] @testset "setindex! and getindex" begin p1 = pts[1] @@ -53,6 +56,7 @@ using RecursiveArrayTools test_vee_hat_from_identity=true, test_inv_diff=true, test_adjoint_inv_diff=true, + test_adjoint_action=true, ) @test isapprox( G, @@ -94,7 +98,15 @@ using RecursiveArrayTools @test compose(G, pts[1], Identity(G)) == pts[1] @test compose(G, Identity(G), pts[1]) == pts[1] - test_group(G, pts, X_pts, X_pts; test_diff=true, test_mutating=false) + test_group( + G, + pts, + X_pts, + X_pts; + test_diff=true, + test_mutating=false, + test_adjoint_action=true, + ) test_manifold(G, pts; is_mutating=false) @test isapprox( G, diff --git a/test/groups/rotation_action.jl b/test/groups/rotation_action.jl index 2a29167258..81db444cef 100644 --- a/test/groups/rotation_action.jl +++ b/test/groups/rotation_action.jl @@ -53,6 +53,8 @@ include("group_utils.jl") atol=atol, ) + @test apply_diff(A_left, Identity(G), m_pts[1], X_pts[1]) === X_pts[1] + @testset "apply_diff_group" begin @test apply_diff_group(A_left, Identity(G), a_X_pts[1], m_pts[1]) ≈ a_X_pts[1] * m_pts[1] diff --git a/test/groups/semidirect_product_group.jl b/test/groups/semidirect_product_group.jl index 83611ca9bd..5b506ca319 100644 --- a/test/groups/semidirect_product_group.jl +++ b/test/groups/semidirect_product_group.jl @@ -4,10 +4,11 @@ include("group_utils.jl") @testset "Semidirect product group" begin M1 = TranslationGroup(2) A = TranslationAction(M1, M1) - G = SemidirectProductGroup(M1, M1, A) + G = SemidirectProductGroup(M1, M1, A, Manifolds.LeftInvariantRepresentation()) @test G === GroupManifold( TranslationGroup(2) × TranslationGroup(2), Manifolds.SemidirectProductOperation(A), + Manifolds.LeftInvariantRepresentation(), ) @test repr(G) == "SemidirectProductGroup($(M1), $(M1), $(A))" @test repr(G.op) == "SemidirectProductOperation($(A))" @@ -28,7 +29,7 @@ include("group_utils.jl") @test p2[G, 1] == p1[M, 1] end - X = log(G, pts[1], pts[1]) + X = log(base_manifold(G), pts[1], pts[1]) Y = zero_vector(G, pts[1]) Z = Manifolds.allocate_result(G, zero_vector, pts[1]) Z = zero_vector!(M, Z, pts[1]) @@ -51,7 +52,7 @@ include("group_utils.jl") eA = identity_element(G) @test isapprox(G, eA, e) @test isapprox(G, e, eA) - W = log(G, eA, pts[1]) - Z = log(G, eA, pts[1]) + W = log(base_manifold(G), eA, pts[1]) + Z = log(base_manifold(G), eA, pts[1]) @test isapprox(G, e, W, Z) end diff --git a/test/groups/special_euclidean.jl b/test/groups/special_euclidean.jl index 949f8752db..3942ccb6e8 100644 --- a/test/groups/special_euclidean.jl +++ b/test/groups/special_euclidean.jl @@ -9,19 +9,33 @@ using Manifolds: LeftForwardAction, LeftBackwardAction, RightForwardAction, RightBackwardAction @testset "Special Euclidean group" begin - for se_parameter in [:field, :type] - @testset "SpecialEuclidean($n)" for n in (2, 3, 4) - G = SpecialEuclidean(n; parameter=se_parameter) + @test repr(SpecialEuclidean(3; vectors=HybridTangentRepresentation())) == + "SpecialEuclidean(3; vectors=HybridTangentRepresentation())" + for (se_parameter, se_vectors) in [ + (:field, LeftInvariantRepresentation()), + (:type, LeftInvariantRepresentation()), + (:field, HybridTangentRepresentation()), + ] + @testset "SpecialEuclidean($n; parameter=$se_parameter, vectors=$se_vectors)" for n in + ( + 2, + 3, + 4, + ) + G = SpecialEuclidean(n; parameter=se_parameter, vectors=se_vectors) if se_parameter === :field @test isa(G, SpecialEuclidean{Tuple{Int}}) else @test isa(G, SpecialEuclidean{TypeParameter{Tuple{n}}}) end - if se_parameter === :field + if se_parameter === :field && se_vectors === LeftInvariantRepresentation() @test repr(G) == "SpecialEuclidean($n; parameter=:field)" - else + elseif se_parameter === :type && se_vectors === LeftInvariantRepresentation() @test repr(G) == "SpecialEuclidean($n)" + elseif se_parameter === :field && se_vectors === HybridTangentRepresentation() + @test repr(G) == + "SpecialEuclidean($n; parameter=:field, vectors=HybridTangentRepresentation())" end M = base_manifold(G) @test M === @@ -92,13 +106,11 @@ using Manifolds: @test affine_matrix(G, Identity(G)) == SDiagonal{n,Float64}(I) w = translate_diff(G, pts[1], Identity(G), X_pts[1]) - w2 = allocate(w) - submanifold_component(w2, 1) .= submanifold_component(w, 1) - submanifold_component(w2, 2) .= - submanifold_component(pts[1], 2) * submanifold_component(w, 2) - w2mat = screw_matrix(G, w2) - @test w2mat ≈ affine_matrix(G, pts[1]) * screw_matrix(G, X_pts[1]) - @test screw_matrix(G, w2mat) === w2mat + if se_vectors isa Manifolds.LeftInvariantRepresentation + w2mat = screw_matrix(G, w) + @test w2mat ≈ screw_matrix(G, X_pts[1]) + @test screw_matrix(G, w2mat) === w2mat + end @test is_vector(G, Identity(G), rand(G; vector_at=Identity(G))) @@ -121,6 +133,7 @@ using Manifolds: basis_types_vecs=basis_types, basis_types_to_from=basis_types, is_mutating=true, + is_tangent_atol_multiplier=1, #test_inplace=true, test_vee_hat=true, exp_log_atol_multiplier=50, @@ -146,6 +159,7 @@ using Manifolds: pts; is_mutating=true, exp_log_atol_multiplier=50, + is_tangent_atol_multiplier=1, test_inner=false, test_norm=false, ) @@ -174,6 +188,7 @@ using Manifolds: basis_types_to_from=basis_types, is_mutating=true, exp_log_atol_multiplier=50, + is_tangent_atol_multiplier=1, ) end end @@ -195,7 +210,7 @@ using Manifolds: pts, X_pts, X_pts; - test_diff=true, + test_diff=true, # fails sometimes test_lie_bracket=true, diff_convs=[(), (LeftForwardAction(),), (RightBackwardAction(),)], atol=1e-9, @@ -206,6 +221,7 @@ using Manifolds: is_mutating=true, #test_inplace=true, test_vee_hat=true, + test_is_tangent=true, # fails exp_log_atol_multiplier=50, ) # specific affine tests @@ -278,66 +294,67 @@ using Manifolds: G = SpecialEuclidean(11) @test affine_matrix(G, Identity(G)) isa Diagonal{Float64,Vector{Float64}} @test affine_matrix(G, Identity(G)) == Diagonal(ones(11)) + end + end - @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 = [ - ArrayPartition(ti, exp(Rn, p, hat(Rn, p, ωi))) for (ti, ωi) in zip(t, ω) - ] - X = ArrayPartition([-1.0, 2.0, 1.0], hat(Rn, p, [1.0, 0.5, -0.5])) - q = ArrayPartition([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 [LeftForwardAction(), RightBackwardAction()] - 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], X_gl, conv), - ) - end + for se_vectors in [LeftInvariantRepresentation(), HybridTangentRepresentation()] + @testset "Explicit embedding in GL(n+1)" begin + G = SpecialEuclidean(3; vectors=se_vectors) + 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 = [ArrayPartition(ti, exp(Rn, p, hat(Rn, p, ωi))) for (ti, ωi) in zip(t, ω)] + X = ArrayPartition([-1.0, 2.0, 1.0], hat(Rn, p, [1.0, 0.5, -0.5])) + q = ArrayPartition([0.0, 0.0, 0.0], p) + + GL = GeneralLinear(4) + SEGL = EmbeddedManifold(G, GL) + @test Manifolds.SpecialEuclideanInGeneralLinear(3; se_vectors=se_vectors) === + 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 [LeftForwardAction(), RightBackwardAction()] + 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], X_gl, conv), + ) end end end @testset "Adjoint action on 𝔰𝔢(3)" begin - G = SpecialEuclidean(3; parameter=:type) + G = SpecialEuclidean(3; parameter=:type, vectors=HybridTangentRepresentation()) 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) @@ -353,9 +370,61 @@ using Manifolds: @test isapprox(G, pts[1], hat(G, pts[1], fXp.data), fXp2) end + @testset "Invariant exp and log, inv_diff" begin + G = SpecialEuclidean(3) + p = ArrayPartition( + [-0.3879800256554809, -1.480242310944754, 0.6859001130634623], + [ + 0.07740722383491305 0.37616397751531383 0.9233140222687134 + -0.8543895289705175 0.5023190688041862 -0.13301911855531234 + -0.513835240601215 -0.7785731918937062 0.36027368816045613 + ], + ) + X = ArrayPartition( + [-0.13357857168916804, -0.6285892085394564, -0.5876201702680527], + [ + 0.0 -0.9387500040623927 -0.02175499382616511 + 0.9387500040623927 0.0 -0.31110293886623697 + 0.02175499382616511 0.31110293886623697 0.0 + ], + ) + q_ref = ArrayPartition( + [-1.188881359581349, -1.7593716755284188, 0.7643154032226447], + [ + 0.4842416235264354 0.37906436972889457 0.788555802493722 + -0.13102724697252144 0.9225290425378667 -0.3630041683300185 + -0.8650675757391926 0.07245943193403166 0.49639472209977575 + ], + ) + q = exp_inv(G, p, X) + @test isapprox(G, q, q_ref) + Y = log_inv(G, p, q) + @test isapprox(G, p, X, Y) + + q2 = similar(q) + exp_inv!(G, q2, p, X) + @test isapprox(G, q, q2) + + Y2 = similar(Y) + log_inv!(G, Y2, p, q) + @test isapprox(G, p, Y, Y2) + + X_inv_diff_ref = ArrayPartition( + [0.8029663714810721, -0.5577248382341342, -1.8086993509926863], + [ + -0.0 0.1952891277629885 0.40160273972404353 + -0.1952891277629885 -5.204170427930421e-17 0.8826592119716654 + -0.4016027397240436 -0.8826592119716655 2.7755575615628914e-17 + ], + ) + @test isapprox(inv_diff(G, p, X), X_inv_diff_ref) + Y3 = similar(Y) + inv_diff!(G, Y3, p, X) + @test isapprox(Y3, X_inv_diff_ref) + end @testset "performance of selected operations" begin for n in [2, 3] - SEn = SpecialEuclidean(n) + SEn = SpecialEuclidean(n; vectors=HybridTangentRepresentation()) Rn = Rotations(n) p = SMatrix{n,n}(I) @@ -382,6 +451,7 @@ using Manifolds: ] end exp(SEn, pts[1], Xs[1]) + exp(base_manifold(SEn), pts[1], Xs[1]) compose(SEn, pts[1], pts[2]) log(SEn, pts[1], pts[2]) log(SEn, pts[1], pts[3]) @@ -394,7 +464,7 @@ using Manifolds: get_vector(SEn, pts[1], csen, DefaultOrthogonalBasis()) # @btime shows 0 but `@allocations` is inaccurate @static if VERSION >= v"1.9-DEV" - @test (@allocations exp(SEn, pts[1], Xs[1])) <= 4 + @test (@allocations exp(base_manifold(SEn), pts[1], Xs[1])) <= 4 @test (@allocations compose(SEn, pts[1], pts[2])) <= 4 if VERSION < v"1.11-DEV" @test (@allocations log(SEn, pts[1], pts[2])) <= 28 diff --git a/test/groups/special_linear.jl b/test/groups/special_linear.jl index c16f508243..80ae54095c 100644 --- a/test/groups/special_linear.jl +++ b/test/groups/special_linear.jl @@ -2,6 +2,8 @@ include("../header.jl") include("group_utils.jl") using NLsolve +using Manifolds: LeftForwardAction, RightBackwardAction + @testset "Special Linear group" begin @testset "basic properties" begin G = SpecialLinear(3) @@ -60,9 +62,8 @@ using NLsolve ) types = [Matrix{Float64}] - pts = - [[2 -1 -3; 4 -1 -6; -1 1 2], [0 2 1; 0 -3 -1; 1 0 2], [-2 0 -1; 1 0 0; -1 -1 2]] - vpts = [[0 -1 -5; 1 2 0; 1 2 -2], [0 -2 1; -2 1 2; -4 2 -1]] + vpts = [[0 -1 -5; 1 2 0; 1 2 -2], [0 -2 1; -2 1 2; -4 2 -1], [0 1 0; 0 0 0; 0 0 0]] + pts = exp.(0.1 .* vpts) retraction_methods = [ Manifolds.GroupExponentialRetraction(LeftForwardAction()), diff --git a/test/header.jl b/test/header.jl index bbffe9d9fc..8be4f238d9 100644 --- a/test/header.jl +++ b/test/header.jl @@ -17,6 +17,7 @@ using DoubleFloats using Quaternions using Random using StaticArrays +using RecursiveArrayTools using Statistics using StatsBase using Test diff --git a/test/manifolds/power_manifold.jl b/test/manifolds/power_manifold.jl index 192df1dced..2e30aee847 100644 --- a/test/manifolds/power_manifold.jl +++ b/test/manifolds/power_manifold.jl @@ -67,54 +67,86 @@ end inverse_retraction_methods = [ManifoldsBase.LogarithmicInverseRetraction()] sphere_dist = Manifolds.uniform_distribution(Ms, @SVector [1.0, 0.0, 0.0]) - power_s1_pt_dist = - Manifolds.PowerPointDistribution(Ms1, sphere_dist, randn(Float64, 3, 5)) - power_s2_pt_dist = - Manifolds.PowerPointDistribution(Ms2, sphere_dist, randn(Float64, 3, 5, 7)) - sphere_tv_dist = - Manifolds.normal_tvector_distribution(Ms, (@MVector [1.0, 0.0, 0.0]), 1.0) - power_s1_tv_dist = Manifolds.PowerFVectorDistribution( - TangentSpace(Ms1, rand(power_s1_pt_dist)), - sphere_tv_dist, - ) - power_s2_tv_dist = Manifolds.PowerFVectorDistribution( - TangentSpace(Ms2, rand(power_s2_pt_dist)), - sphere_tv_dist, - ) - - id_rot = @SMatrix [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0] - rotations_dist = Manifolds.normal_rotation_distribution(Mr, id_rot, 1.0) - power_r1_pt_dist = - Manifolds.PowerPointDistribution(Mr1, rotations_dist, randn(Float64, 3, 3, 5)) - power_rn1_pt_dist = Manifolds.PowerPointDistribution( - Mrn1, - rotations_dist, - [randn(Float64, 3, 3) for i in 1:5], - ) - power_r2_pt_dist = - Manifolds.PowerPointDistribution(Mr2, rotations_dist, randn(Float64, 3, 3, 5, 7)) - power_rn2_pt_dist = Manifolds.PowerPointDistribution( - Mrn2, - rotations_dist, - [randn(Float64, 3, 3) for i in 1:5, j in 1:7], - ) - rotations_tv_dist = Manifolds.normal_tvector_distribution(Mr, MMatrix(id_rot), 1.0) - power_r1_tv_dist = Manifolds.PowerFVectorDistribution( - TangentSpace(Mr1, rand(power_r1_pt_dist)), - rotations_tv_dist, - ) - power_rn1_tv_dist = Manifolds.PowerFVectorDistribution( - TangentSpace(Mrn1, rand(power_rn1_pt_dist)), - rotations_tv_dist, - ) - power_r2_tv_dist = Manifolds.PowerFVectorDistribution( - TangentSpace(Mr2, rand(power_r2_pt_dist)), - rotations_tv_dist, - ) - power_rn2_tv_dist = Manifolds.PowerFVectorDistribution( - TangentSpace(Mrn2, rand(power_rn2_pt_dist)), - rotations_tv_dist, - ) + + point_distributions_Mrn1 = [] + tvector_distributions_Mrn1 = [] + point_distributions_Mr2 = [] + tvector_distributions_Mr2 = [] + point_distributions_Mrn2 = [] + tvector_distributions_Mrn2 = [] + point_distributions_Ms2 = [] + tvector_distributions_Ms2 = [] + if VERSION >= v"1.9" + PowerPointDistribution = + Base.get_extension(Manifolds, :ManifoldsDistributionsExt).PowerPointDistribution + + PowerFVectorDistribution = + Base.get_extension( + Manifolds, + :ManifoldsDistributionsExt, + ).PowerFVectorDistribution + power_s1_pt_dist = PowerPointDistribution(Ms1, sphere_dist, randn(Float64, 3, 5)) + power_s2_pt_dist = PowerPointDistribution(Ms2, sphere_dist, randn(Float64, 3, 5, 7)) + push!(point_distributions_Ms2, power_s2_pt_dist) + + sphere_tv_dist = + Manifolds.normal_tvector_distribution(Ms, (@MVector [1.0, 0.0, 0.0]), 1.0) + power_s1_tv_dist = PowerFVectorDistribution( + TangentSpace(Ms1, rand(power_s1_pt_dist)), + sphere_tv_dist, + ) + power_s2_tv_dist = PowerFVectorDistribution( + TangentSpace(Ms2, rand(power_s2_pt_dist)), + sphere_tv_dist, + ) + push!(tvector_distributions_Ms2, power_s2_tv_dist) + + id_rot = @SMatrix [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0] + rotations_dist = Manifolds.normal_rotation_distribution(Mr, id_rot, 1.0) + power_r1_pt_dist = + PowerPointDistribution(Mr1, rotations_dist, randn(Float64, 3, 3, 5)) + power_rn1_pt_dist = PowerPointDistribution( + Mrn1, + rotations_dist, + [randn(Float64, 3, 3) for i in 1:5], + ) + push!(point_distributions_Mrn1, power_rn1_pt_dist) + power_r2_pt_dist = + PowerPointDistribution(Mr2, rotations_dist, randn(Float64, 3, 3, 5, 7)) + power_rn2_pt_dist = PowerPointDistribution( + Mrn2, + rotations_dist, + [randn(Float64, 3, 3) for i in 1:5, j in 1:7], + ) + push!(point_distributions_Mr2, power_r2_pt_dist) + rotations_tv_dist = Manifolds.normal_tvector_distribution(Mr, MMatrix(id_rot), 1.0) + power_r1_tv_dist = PowerFVectorDistribution( + TangentSpace(Mr1, rand(power_r1_pt_dist)), + rotations_tv_dist, + ) + power_rn1_tv_dist = PowerFVectorDistribution( + TangentSpace(Mrn1, rand(power_rn1_pt_dist)), + rotations_tv_dist, + ) + push!(tvector_distributions_Mrn1, power_rn1_tv_dist) + power_r2_tv_dist = PowerFVectorDistribution( + TangentSpace(Mr2, rand(power_r2_pt_dist)), + rotations_tv_dist, + ) + power_rn2_tv_dist = PowerFVectorDistribution( + TangentSpace(Mrn2, rand(power_rn2_pt_dist)), + rotations_tv_dist, + ) + push!(tvector_distributions_Mrn2, power_rn2_tv_dist) + + MPointSupport = + Base.get_extension(Manifolds, :ManifoldsDistributionsExt).MPointSupport + FVectorSupport = + Base.get_extension(Manifolds, :ManifoldsDistributionsExt).FVectorSupport + + Test.@test Distributions.support(power_s1_pt_dist) isa MPointSupport + Test.@test Distributions.support(power_s1_tv_dist) isa FVectorSupport + end @testset "get_component, set_component!, getindex and setindex!" begin p1 = randn(3, 5) @@ -172,7 +204,7 @@ end basis_types = (DefaultOrthonormalBasis(), ProjectedOrthonormalBasis(:svd)) for T in types_s2 @testset "Type $(trim(string(T)))..." begin - pts2 = [convert(T, rand(power_s2_pt_dist)) for _ in 1:3] + pts2 = [convert(T, rand(Ms2)) for _ in 1:3] test_manifold( Ms2, pts2; @@ -181,8 +213,8 @@ end test_vee_hat=true, retraction_methods=retraction_methods, inverse_retraction_methods=inverse_retraction_methods, - point_distributions=[power_s2_pt_dist], - tvector_distributions=[power_s2_tv_dist], + point_distributions=point_distributions_Ms2, + tvector_distributions=tvector_distributions_Ms2, rand_tvector_atol_multiplier=6.0, retraction_atol_multiplier=12, is_tangent_atol_multiplier=12.0, @@ -194,7 +226,7 @@ end for T in types_rn1 @testset "Type $(trim(string(T)))..." begin - pts1 = [convert(T, rand(power_rn1_pt_dist)) for _ in 1:3] + pts1 = [convert(T, rand(Mrn1)) for _ in 1:3] test_manifold( Mrn1, pts1; @@ -203,8 +235,8 @@ end test_vee_hat=true, retraction_methods=retraction_methods, inverse_retraction_methods=inverse_retraction_methods, - point_distributions=[power_rn1_pt_dist], - tvector_distributions=[power_rn1_tv_dist], + point_distributions=point_distributions_Mrn1, + tvector_distributions=tvector_distributions_Mrn1, basis_types_to_from=basis_types, rand_tvector_atol_multiplier=500.0, retraction_atol_multiplier=12, @@ -218,7 +250,7 @@ end end for T in types_r2 @testset "Type $(trim(string(T)))..." begin - pts2 = [convert(T, rand(power_r2_pt_dist)) for _ in 1:3] + pts2 = [convert(T, rand(Mr2)) for _ in 1:3] test_manifold( Mr2, pts2; @@ -227,8 +259,8 @@ end test_vee_hat=true, retraction_methods=retraction_methods, inverse_retraction_methods=inverse_retraction_methods, - point_distributions=[power_r2_pt_dist], - tvector_distributions=[power_r2_tv_dist], + point_distributions=point_distributions_Mr2, + tvector_distributions=tvector_distributions_Mr2, rand_tvector_atol_multiplier=8.0, retraction_atol_multiplier=12, is_tangent_atol_multiplier=12.0, @@ -239,7 +271,7 @@ end end for T in types_rn2 @testset "Type $(trim(string(T)))..." begin - pts2 = [convert(T, rand(power_rn2_pt_dist)) for _ in 1:3] + pts2 = [convert(T, rand(Mrn2)) for _ in 1:3] test_manifold( Mrn2, pts2; @@ -248,8 +280,8 @@ end test_vee_hat=true, retraction_methods=retraction_methods, inverse_retraction_methods=inverse_retraction_methods, - point_distributions=[power_rn2_pt_dist], - tvector_distributions=[power_rn2_tv_dist], + point_distributions=point_distributions_Mrn2, + tvector_distributions=tvector_distributions_Mrn2, rand_tvector_atol_multiplier=8.0, retraction_atol_multiplier=12, is_tangent_atol_multiplier=12.0, diff --git a/test/manifolds/product_manifold.jl b/test/manifolds/product_manifold.jl index 9f12478ea0..eb548a2ee1 100644 --- a/test/manifolds/product_manifold.jl +++ b/test/manifolds/product_manifold.jl @@ -120,8 +120,6 @@ using RecursiveArrayTools: ArrayPartition M2, Distributions.MvNormal(zero(pts_r2[1]), 1.0 * I), ) - distr_tv_M1 = Manifolds.normal_tvector_distribution(M1, pts_sphere[1], 1.0) - distr_tv_M2 = Manifolds.normal_tvector_distribution(M2, pts_r2[1], 1.0) @test injectivity_radius(Mse, pts[1]) ≈ π @test injectivity_radius(Mse) ≈ π @test injectivity_radius(Mse, pts[1], ExponentialRetraction()) ≈ π @@ -143,13 +141,45 @@ using RecursiveArrayTools: ArrayPartition inverse_retract(Mse, pts[1], pts[2], default_inverse_retraction_method(Mse)), ) + Mse_point_distributions = [] + Mse_tvector_distributions = [] + + if VERSION >= v"1.9" + distr_tv_M1 = Manifolds.normal_tvector_distribution(M1, pts_sphere[1], 1.0) + distr_tv_M2 = Manifolds.normal_tvector_distribution(M2, pts_r2[1], 1.0) + + ProductPointDistribution = + Base.get_extension( + Manifolds, + :ManifoldsDistributionsExt, + ).ProductPointDistribution + push!(Mse_point_distributions, ProductPointDistribution(distr_M1, distr_M2)) + + ProductFVectorDistribution = + Base.get_extension( + Manifolds, + :ManifoldsDistributionsExt, + ).ProductFVectorDistribution + push!( + Mse_tvector_distributions, + ProductFVectorDistribution(distr_tv_M1, distr_tv_M2), + ) + + MPointSupport = + Base.get_extension(Manifolds, :ManifoldsDistributionsExt).MPointSupport + FVectorSupport = + Base.get_extension(Manifolds, :ManifoldsDistributionsExt).FVectorSupport + + Test.@test Distributions.support(Mse_point_distributions[1]) isa MPointSupport + Test.@test Distributions.support(Mse_tvector_distributions[1]) isa + FVectorSupport + end + test_manifold( Mse, pts; - point_distributions=[Manifolds.ProductPointDistribution(distr_M1, distr_M2)], - tvector_distributions=[ - Manifolds.ProductFVectorDistribution(distr_tv_M1, distr_tv_M2), - ], + point_distributions=Mse_point_distributions, + tvector_distributions=Mse_tvector_distributions, test_injectivity_radius=true, test_musical_isomorphisms=true, musical_isomorphism_bases=[DefaultOrthonormalBasis()], @@ -193,7 +223,7 @@ using RecursiveArrayTools: ArrayPartition end @testset "empty allocation" begin - p = allocate_result(Mse, uniform_distribution) + p = allocate_result(Mse, Manifolds.uniform_distribution) @test isa(p, ArrayPartition) @test size(p[Mse, 1]) == (3,) @test size(p[Mse, 2]) == (2,) @@ -201,9 +231,9 @@ using RecursiveArrayTools: ArrayPartition @testset "Uniform distribution" begin Mss = ProductManifold(Sphere(2), Sphere(2)) - p = rand(uniform_distribution(Mss)) + p = rand(Manifolds.uniform_distribution(Mss)) @test is_point(Mss, p) - @test is_point(Mss, rand(uniform_distribution(Mss, p))) + @test is_point(Mss, rand(Manifolds.uniform_distribution(Mss, p))) end @testset "Atlas & Induced Basis" begin diff --git a/test/manifolds/rotations.jl b/test/manifolds/rotations.jl index a36b646597..a19ac5e951 100644 --- a/test/manifolds/rotations.jl +++ b/test/manifolds/rotations.jl @@ -102,6 +102,8 @@ include("../header.jl") gtsd_mvector = Manifolds.normal_tvector_distribution(M, (@MMatrix [1.0 0.0; 0.0 1.0]), 1.0) @test isa(rand(gtsd_mvector), MMatrix) + + Test.@test Distributions.support(usd_mmatrix).manifold == M end Random.seed!(42) diff --git a/test/manifolds/sphere.jl b/test/manifolds/sphere.jl index 643f99c968..aaa4f39ad9 100644 --- a/test/manifolds/sphere.jl +++ b/test/manifolds/sphere.jl @@ -248,7 +248,7 @@ using ManifoldsBase: TFVector @testset "Distributions" begin sphere = Sphere(2) - haar_measure = uniform_distribution(sphere) + haar_measure = Manifolds.uniform_distribution(sphere) pts = rand(haar_measure, 5) @test all(p -> is_point(sphere, p), pts) end diff --git a/test/manifolds/symplecticgrassmann.jl b/test/manifolds/symplecticgrassmann.jl index 43accdf8df..b2b77c709f 100644 --- a/test/manifolds/symplecticgrassmann.jl +++ b/test/manifolds/symplecticgrassmann.jl @@ -36,8 +36,8 @@ include("../header.jl") -0.02823014 0.00029946 -0.04196034 -0.04145413 ] @testset "Basics" begin - @test repr(M) == "SymplecticGrassmann(6, 4; field=ℝ)" - @test repr(Mf) == "SymplecticGrassmann(6, 4; field=ℝ; parameter=:field)" + @test repr(M) == "SymplecticGrassmann(6, 4)" + @test repr(Mf) == "SymplecticGrassmann(6, 4; parameter=:field)" @test manifold_dimension(M) == 4 * (6 - 4) for _M in [M, Mf] @test is_point(M, p) diff --git a/test/manifolds/symplecticstiefel.jl b/test/manifolds/symplecticstiefel.jl index bee161cd0b..23250d153e 100644 --- a/test/manifolds/symplecticstiefel.jl +++ b/test/manifolds/symplecticstiefel.jl @@ -138,7 +138,7 @@ end @testset "Basics" begin @test_throws ArgumentError SymplecticStiefel(5, 4) @test_throws ArgumentError SymplecticStiefel(6, 3) - @test repr(M) == "SymplecticStiefel(6, 4; field=ℝ)" + @test repr(M) == "SymplecticStiefel(6, 4)" @test representation_size(M) == (6, 4) @test base_manifold(M) === M @test get_total_space(M) == SymplecticMatrices(6) @@ -306,7 +306,7 @@ end @testset "field parameter" begin M = SymplecticStiefel(6, 4; parameter=:field) @test typeof(get_embedding(M)) === Euclidean{Tuple{Int,Int},ℝ} - @test repr(M) == "SymplecticStiefel(6, 4; field=ℝ; parameter=:field)" + @test repr(M) == "SymplecticStiefel(6, 4; parameter=:field)" @test get_total_space(M) == SymplecticMatrices(6; parameter=:field) end end diff --git a/test/metric.jl b/test/metric.jl index 6441406310..c19fbb2b82 100644 --- a/test/metric.jl +++ b/test/metric.jl @@ -192,12 +192,19 @@ function Manifolds.get_vector_orthonormal!( ) return Y .= X end + Manifolds.is_default_metric(::BaseManifold, ::DefaultBaseManifoldMetric) = true -function Manifolds.projected_distribution(M::BaseManifold, d) - return ProjectedPointDistribution(M, d, project!, rand(d)) -end -function Manifolds.projected_distribution(M::BaseManifold, d, p) - return ProjectedPointDistribution(M, d, project!, p) + +if VERSION >= v"1.9" + const ProjectedPointDistribution = + Base.get_extension(Manifolds, :ManifoldsDistributionsExt).ProjectedPointDistribution + + function Manifolds.projected_distribution(M::BaseManifold, d) + return ProjectedPointDistribution(M, d, project!, rand(d)) + end + function Manifolds.projected_distribution(M::BaseManifold, d, p) + return ProjectedPointDistribution(M, d, project!, p) + end end function Manifolds.flat!( @@ -555,8 +562,10 @@ Manifolds.inner(::MetricManifold{ℝ,<:AbstractManifold{ℝ},Issue539Metric}, p, @test_throws MethodError local_metric_jacobian(MM2, p, B_p) @test_throws MethodError christoffel_symbols_second_jacobian(MM2, p, B_p) # MM falls back to nondefault error - @test_throws MethodError projected_distribution(MM, 1, p) - @test_throws MethodError projected_distribution(MM, 1) + if VERSION >= v"1.9" + @test_throws MethodError Manifolds.projected_distribution(MM, 1, p) + @test_throws MethodError Manifolds.projected_distribution(MM, 1) + end @test inner(MM2, p, X, Y) === inner(M, p, X, Y) @test norm(MM2, p, X) === norm(M, p, X) @@ -595,16 +604,18 @@ Manifolds.inner(::MetricManifold{ℝ,<:AbstractManifold{ℝ},Issue539Metric}, p, @test is_point(MM2, p) === is_point(M, p) @test is_vector(MM2, p, X) === is_vector(M, p, X) - a = Manifolds.projected_distribution( - M, - Distributions.MvNormal(zero(zeros(3)), 1.0 * I), - ) - b = Manifolds.projected_distribution( - MM2, - Distributions.MvNormal(zero(zeros(3)), 1.0 * I), - ) - @test isapprox(Matrix(a.distribution.Σ), Matrix(b.distribution.Σ)) - @test isapprox(a.distribution.μ, b.distribution.μ) + if VERSION >= v"1.9" + a = Manifolds.projected_distribution( + M, + Distributions.MvNormal(zero(zeros(3)), 1.0 * I), + ) + b = Manifolds.projected_distribution( + MM2, + Distributions.MvNormal(zero(zeros(3)), 1.0 * I), + ) + @test isapprox(Matrix(a.distribution.Σ), Matrix(b.distribution.Σ)) + @test isapprox(a.distribution.μ, b.distribution.μ) + end @test get_basis(M, p, DefaultOrthonormalBasis()).data == get_basis(MM2, p, DefaultOrthonormalBasis()).data @test_throws MethodError get_basis(MM, p, DefaultOrthonormalBasis()) diff --git a/test/statistics.jl b/test/statistics.jl index 7cbb10b753..d83d234254 100644 --- a/test/statistics.jl +++ b/test/statistics.jl @@ -756,7 +756,7 @@ end @test m == mf μ = project(R, randn(3, 3)) - d = normal_tvector_distribution(R, μ, 0.1) + d = Manifolds.normal_tvector_distribution(R, μ, 0.1) x = [exp(R, μ, rand(rng, d)) for _ in 1:10] w = pweights([rand(rng) for _ in 1:length(x)]) m = mean(R, x, w) diff --git a/test/test_deprecated.jl b/test/test_deprecated.jl index 926f9c0720..3e74a0d5b8 100644 --- a/test/test_deprecated.jl +++ b/test/test_deprecated.jl @@ -1,33 +1,3 @@ using Manifolds, ManifoldsBase, Random, Test -using StatsBase: AbstractWeights, pweights -using Random: GLOBAL_RNG, seed! - -@testset "Deprecation tests" begin - @testset "Deprecate extrinsic_method= keyword" begin - rng = MersenneTwister(47) - S = Sphere(2) - x = [normalize(randn(rng, 3)) for _ in 1:10] - w = pweights([rand(rng) for _ in 1:length(x)]) - mg1 = mean(S, x, w, ExtrinsicEstimation(EfficientEstimator())) - # Statistics 414-418, deprecate former extrinsic_method keyword - mg2 = mean( - S, - x, - w, - ExtrinsicEstimation(EfficientEstimator()); - extrinsic_method=EfficientEstimator(), - ) - @test isapprox(S, mg1, mg2) - mg3 = median(S, x, w, ExtrinsicEstimation(CyclicProximalPointEstimation())) - # Statistics 692-696, deprecate former extrinsic_method keyword - mg4 = median( - S, - x, - w, - ExtrinsicEstimation(CyclicProximalPointEstimation()); - extrinsic_method=CyclicProximalPointEstimation(), - ) - @test isapprox(S, mg3, mg4) - end -end +@testset "Deprecation tests" begin end diff --git a/tutorials/Project.toml b/tutorials/Project.toml index 3d0f29a1a0..98bf338cff 100644 --- a/tutorials/Project.toml +++ b/tutorials/Project.toml @@ -4,6 +4,7 @@ BoundaryValueDiffEq = "764a87c0-6b3e-53db-9096-fe964310641d" CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def" +HybridArrays = "1baab800-613f-4b0a-84e4-9cd3431bfbb9" IJulia = "7073ff75-c697-5162-941a-fcdaad2a7d2a" Manifolds = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e" ManifoldsBase = "3362f125-f0bb-47a3-aa74-596ffd7ef2fb" @@ -20,8 +21,10 @@ BoundaryValueDiffEq = "4, 5" CSV = "0.10" DataFrames = "1" Distributions = "0.22.6, 0.23, 0.24, 0.25" +HybridArrays = "0.4" IJulia = "1" -Manifolds = "0.9.0" +Manifolds = "0.10.0" ManifoldsBase = "0.15.0" MultivariateStats = "0.10" +RecursiveArrayTools = "3" StaticArrays = "1" diff --git a/tutorials/_quarto.yml b/tutorials/_quarto.yml index 5befbaf19d..0966a12e28 100644 --- a/tutorials/_quarto.yml +++ b/tutorials/_quarto.yml @@ -21,7 +21,7 @@ execute: format: commonmark: - variant: -raw_html+tex_math_dollars + variant: -raw_html+tex_math_dollars+pipe_tables+footnotes wrap: preserve jupyter: julia-1.10 diff --git a/tutorials/getstarted.qmd b/tutorials/getstarted.qmd index a0d8d3d7d7..1194c16bbe 100644 --- a/tutorials/getstarted.qmd +++ b/tutorials/getstarted.qmd @@ -224,7 +224,8 @@ M₆ = M₂ × M₃ Since now the representations might differ from element to element, we have to encapsulate these in their own type. ```{julia} -p₃ = Manifolds.ArrayPartition([0, 0, 1], [0, 1, 0]) +using RecursiveArrayTools: ArrayPartition +p₃ = ArrayPartition([0, 0, 1], [0, 1, 0]) ``` Here `ArrayPartition` taken from [🔗 `RecursiveArrayTools.jl`](https://github.com/SciML/RecursiveArrayTools.jl) to store the point on the product manifold efficiently in one array, still allowing efficient access to the product elements. diff --git a/tutorials/groups.qmd b/tutorials/groups.qmd new file mode 100644 index 0000000000..3f2c207fe4 --- /dev/null +++ b/tutorials/groups.qmd @@ -0,0 +1,238 @@ +--- +title: work with groups +toc: true +--- + +```{julia} +#| echo: false +#| code-fold: true +#| output: false +using Pkg; +cd(@__DIR__) +Pkg.activate("."); # for reproducibility use the local tutorial environment. +Pkg.develop(path="../") # a trick to work on the local dev version +using Markdown +``` + +This is a short overview of group support in [`Manifolds.jl`](https://juliamanifolds.github.io/Manifolds.jl/) and how to get started working with them. + +Groups currently available in `Manifolds.jl` are listed in ``[group section](@ref GroupManifoldSection)``{=commonmark}. + +You can read more about the theory of Lie groups for example in [Chirikjian:2012](@cite). +An example application of Lie groups in robotics is described in [SolaDerayAtchuthan:2021](@cite). + +First, let's load libraries we will use. +[`RecursiveArrayTools.jl`](https://github.com/SciML/RecursiveArrayTools.jl) is necessary because its `ArrayPartition` is used as one of the possible representations of elements of product and semidirect product groups. +[`StaticArrays.jl`](https://github.com/JuliaArrays/StaticArrays.jl) can be used to speed up some operations. + +```{julia} +using Manifolds, RecursiveArrayTools, StaticArrays +``` + +## Introduction: group of rotations on a plane + +Let's first consider an example of the group of rotations of a plane, $\operatorname{SO}(2)$. +They can be represented in several ways, for example as angles of rotation (which corresponds to [`RealCircleGroup`](@ref)), unit complex numbers ([`CircleGroup`](@ref)), or rotation matrices ([`SpecialOrthogonal`](@ref)). +Let's consider the last representation since it is the most nontrivial one and can be more easily generalized to other groups. + +The associated manifolds and groups are defined by: + +```{julia} +G = SpecialOrthogonal(2) +M = base_manifold(G) +@assert M === Rotations(2) +``` + +This duality (Lie group and the underlying manifold being separate) is a common pattern in `Manifolds.jl`. The group `G` can be used for both Lie group-specific operations and metric-specific operation, while the manifold `M` only allows using manifold and metric operations. This way groups can be specialized in ways not relevant to plain manifolds, and if someone doesn't use the groups structure, they don't have to consider it by just using the manifold. + +Some basic definitions + +```{julia} +# default basis +B = DefaultOrthogonalBasis() +# Identity rotation +p0 = @SMatrix [1.0 0; 0 1] + +# Group identity element of a special type +IG = Identity(G) +``` + + +Let's say we want to define a manifold point `p_i` some rotation θ from the [`identity_element`](@ref) reference rotation `p0` (another point on the manifold that we will use as reference) + +```{julia} +# + radians rotation from x-axis on plane to point i +xθi = π/6 +``` + +### From Coordinates + +To get our first Lie algebra element we can use the [`hat`](@ref) function which is commonly used in robotics, or equivalently a more generalized [`get_vector`](@ref), function: + +```{julia} +X_ = hat(G, IG, xθi) # specific definition to Lie algebras +xXi = get_vector(G, p0, xθi, B) # generalized definition beyond Lie algebras +println(xXi) +@assert isapprox(X_, xXi) +``` + +Note that `hat` here assumes a default (orthogonal) basis for the more general `get_vector`. + +::: {.callout-note} +In this case, the same would work given the base manifold [`Rotations(2)`](@ref): +```{julia} +_X_ = hat(M, p0, xθi) # Lie groups definition +_X = get_vector(M, p0, xθi, B) # generalized definition +@assert _X_ == xXi; @assert _X == xXi +``` +One more caveat here is that for the Rotation matrices, the tangent vectors are always stored as elements from the Lie algebra. +::: + +Now, we can transform this algebra element to a point on the manifold using the exponential map [`exp`](@ref): +```{julia} +xRi = exp(G, p0, xXi) +# similarly for known underlying manifold +xRi_ = exp(M, p0, xXi) + +@assert isapprox(xRi, xRi_) +``` + +### To Coordinates + +The logarithmic map transforms elements of the group back to its Lie algebra: +```{julia} +xXi_ = log(G, p0, xRi) +xXi__ = log(M, p0, xRi) +@assert xXi ≈ xXi__ +``` + +Similarly, the coordinate values can be extracted from the algebra elements using [`vee`](@ref), or using the more generalized [`get_coordinates`](@ref): + +```{julia} +# extracting coordinates using vee +xθi__ = vee(G, p0, xXi_)[1] +_xθi__ = vee(M, p0, xXi_)[1] + +# OR, the preferred generalized get_coordinate function +xθi_ = get_coordinates(G, p0, xXi_, B)[1] +_xθi_ = get_coordinates(M, p0, xXi_, B)[1] + +# confirm all versions are correct +@assert isapprox(xθi, xθi_); @assert isapprox(xθi, _xθi_) +@assert isapprox(xθi, xθi__); @assert isapprox(xθi, _xθi__) +``` + + +### Actions and Operations + +With the basics in hand on how to move between the coordinate, algebra, and group representations, let's briefly look at composition and application of points on the manifold. For example, a [`Rotations`](@ref) manifold is the mathematical representation, but the points have an application purpose in retaining information regarding a specific rotation. + +Points from a Lie group may have an associated action (i.e. a rotation) which we [`apply`](@ref). Consider rotating through `θ = π/6` three vectors `V` from their native domain `Euclidean(2)`, from the reference point `a` to a new point `b`. Engineering disciplines sometimes refer to the action of a manifold point `a` or `b` as reference frames. More generally, by taking the tangent space at point `p`, we are defining a local coordinate frame with basis `B`, and should not be confused with "reference frame" `a` or `b`. + +Keeping with our two-dimensional example above: +```{julia} +aV1 = [1; 0] +aV2 = [0; 1] +aV3 = [10; 10] + +A_left = RotationAction(Euclidean(2), G) + +bθa = π/6 +bXa = get_vector(base_manifold(G), p0, bθa, B) + +bRa = exp(G, p0, bXa) + +for aV in [aV1; aV2; aV3] + bV = apply(A_left, bRa, aV) + # test we are getting the rotated vectors in Euclidean(2) as expected + @assert isapprox(bV[1], norm(aV) * cos(bθa)) + @assert isapprox(bV[2], norm(aV) * sin(bθa)) +end +``` + +::: {.callout-note} +In general, actions are usually non-commutative and the user must therefore be aware whether they want to use [`LeftAction`](@ref) or [`RightAction`](@ref). In this case, the default `LeftAction()` is used. +::: + +Finally, the actions (i.e. points from a manifold) can be [`compose`](@ref)d together. Consider putting together two rotations `aRb` and `bRc` such that a single composite rotation `aRc` is found. The next bit of code composes five rotations of `π/4` increments: +```{julia} +A_left = RotationAction(M, G) +aRi = copy(p0) + +iθi_ = π/4 +x_θ = get_vector(M, p0, iθi_, B) #hat(Rn, R0, θ) +iRi_ = exp(M, p0, x_θ) + +# do 5 times over: +# Ri_ = Ri*iRi_ +for i in 1:5 + aRi = compose(A_left, aRi, iRi_) +end + +# drop back to a algebra, then coordinate +aXi = log(G, p0, aRi) +aθi = get_coordinates(G, p0, aXi, B) + +# should wrap around to 3rd quadrant of xy-plane +@assert isapprox(-3π/4, aθi[1]) +``` + +::: {.callout-warning} +`compose` or `apply` must be done with group (not algebra) elements. This example shows how these two element types can easily be confused, since both the manifold group and algebra elements can have exactly the same data storage type -- i.e. a 2x2 matrix. +::: + +As a last note, other rotation representations, including quaternions, Pauli matrices, etc., have similar features. A contrasting example in rotations, however, are Euler angles which can also store rotation information but quickly becomes problematic with familiar problems such as ["gimbal-lock"](https://en.wikipedia.org/wiki/Gimbal_lock). + +## Relationship between groups, metrics and connections + +Group structure provides a canonical way to define [📖 exponential](https://en.wikipedia.org/wiki/Exponential_map_(Lie_theory)) and logarithmic maps from the Lie algebra. +They can be calculated in `Manifolds.jl` using the [`exp_lie`](@ref) and [`log_lie`](@ref) functions. Such exponential and logarithmic maps can be extended invariantly to tangent spaces at any point of the Lie group. This extension is implemented using functions [`exp_inv`](@ref) and [`log_inv`](@ref). + +Finally, there are `log` and `exp` functions which are metric (or connection)-related functions in `Manifolds.jl`. For groups which can be equipped with a bi-invariant metric, `log` and `log_inv` return the same result, similarly `exp` and `exp_inv`. However, only compact groups and their products with Euclidean spaces can have a bi-invariant metric (see for example Theorem 21.9 in [GallierQuaintance:2020](@cite)). A prominent example of a Lie group without a bi-invariant metric is the special Euclidean group (in two or more dimensions). Then we have a choice between a metric but non-invariant exponential map (which is generally the default choice for `exp`) or a non-metric, invariant exponential map (`exp_inv`). Which one should be used depends on whether being metric or being invariant is more important in a particular application. + +```{julia} +G = SpecialEuclidean(2) +p = ArrayPartition([1.0, -1.0], xRi) +X = ArrayPartition([2.0, -3.0], aXi) +q_m = exp(G, p, X) +println(q_m) +q_i = exp_inv(G, p, X) +println(q_i) + +``` + +As we can see, the results differ. We can observe the invariance as follows: +```{julia} +p2 = ArrayPartition([2.0, -1.0], xRi) +q1_m = exp(G, translate(G, p2, p), translate_diff(G, p2, p, X)) +q2_m = translate(G, p2, exp(G, p, X)) +println(isapprox(q1_m, q2_m)) + +q1_i = exp_inv(G, translate(G, p2, p), translate_diff(G, p2, p, X)) +q2_i = translate(G, p2, exp_inv(G, p, X)) +println(isapprox(q1_i, q2_i)) +``` + +Now, `q1_m` and `q2_m` are different due to non-invariance of the metric connection but `q1_i` and `q2_i` are equal due to invarianced of `exp_inv`. + +The following table outlines invariance of `exp` and `log` of various groups. + +| Group | Zero torsion connection | Invariant | +|:---|:---|:---:| +| `ProductGroup` | Product of connections in each submanifold | 🟡^[Yes if all component connections are invariant separately, otherwise no] | +| `SemidirectProductGroup` | Same as underlying product | ❌ | +| `TranslationGroup` | `CartanSchoutenZero` | ✅ | +| `CircleGroup` | `CartanSchoutenZero` | ✅ | +| `GeneralLinearGroup` | Metric connection from the left invariant metric induced from the standard basis on the Lie algebra | ❌ | +| `GeneralUnitaryMultiplicationGroup` | `CartanSchoutenZero` (explicitly) | ✅ | +| `HeisenbergGroup` | `CartanSchoutenZero` | ✅ | +| `SpecialLinearGroup` | Same as `GeneralLinear` | ❌ | + +## Literature + +````{=commonmark} +```@bibliography +Pages = ["groups.md"] +Canonical=false +``` +```` \ No newline at end of file