From 5e7f6ae76e91b96755f282ae202fb84d77fa2614 Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Mon, 16 Oct 2023 10:36:28 +0200 Subject: [PATCH] Define Meta Manifolds already in ManifoldsBase (#169) * Move Tangent space and Vector Bundle to ManifoldsBase. * add Requires.jl * Work on Test coverage. * Switch do Documenter.jl 1.0 and introduce DocumenterCitations. * update CI * fix tests * Some ProductManifold fixes * run ProductManifold tests on CI * tests and fixes * Documentation rework I * update PowerManifold to optionally static * update NEWS --------- Co-authored-by: Mateusz Baran --- .github/workflows/ci-nightly.yml | 25 - .github/workflows/ci.yml | 5 +- NEWS.md | 28 +- Project.toml | 14 +- docs/Project.toml | 6 +- docs/make.jl | 9 +- docs/src/assets/citations.css | 19 + docs/src/manifolds.md | 37 +- docs/src/metamanifolds.md | 37 + docs/src/references.bib | 87 ++ docs/src/references.md | 4 + docs/src/types.md | 7 +- docs/src/vector_transports.md | 2 +- .../ManifoldsBaseRecursiveArrayToolsExt.jl | 69 + .../ProductManifoldRecursiveArrayToolsExt.jl | 375 +++++ src/DefaultManifold.jl | 8 + src/Fiber.jl | 51 + src/ManifoldsBase.jl | 66 +- src/PowerManifold.jl | 189 ++- src/ProductManifold.jl | 1212 +++++++++++++++++ src/TangentSpace.jl | 254 ++++ src/VectorFiber.jl | 11 + src/bases.jl | 45 +- src/point_vector_fallbacks.jl | 20 + src/retractions.jl | 54 +- src/vector_spaces.jl | 7 +- src/vector_transport.jl | 46 +- test/bases.jl | 55 +- test/default_manifold.jl | 11 + test/fibers.jl | 101 ++ test/power.jl | 63 +- test/product_manifold.jl | 560 ++++++++ test/runtests.jl | 30 +- test/test_manifolds.jl | 226 +++ test/test_sphere.jl | 55 +- 35 files changed, 3512 insertions(+), 276 deletions(-) delete mode 100644 .github/workflows/ci-nightly.yml create mode 100644 docs/src/assets/citations.css create mode 100644 docs/src/metamanifolds.md create mode 100644 docs/src/references.bib create mode 100644 docs/src/references.md create mode 100644 ext/ManifoldsBaseRecursiveArrayToolsExt/ManifoldsBaseRecursiveArrayToolsExt.jl create mode 100644 ext/ManifoldsBaseRecursiveArrayToolsExt/ProductManifoldRecursiveArrayToolsExt.jl create mode 100644 src/Fiber.jl create mode 100644 src/ProductManifold.jl create mode 100644 src/TangentSpace.jl create mode 100644 src/VectorFiber.jl create mode 100644 test/fibers.jl create mode 100644 test/product_manifold.jl create mode 100644 test/test_manifolds.jl diff --git a/.github/workflows/ci-nightly.yml b/.github/workflows/ci-nightly.yml deleted file mode 100644 index c68926b4..00000000 --- a/.github/workflows/ci-nightly.yml +++ /dev/null @@ -1,25 +0,0 @@ -name: CI-nightly -on: - push: - branches: [master] - tags: [v*] - pull_request: - -jobs: - test: - name: Julia nightly - ${{ matrix.os }} - runs-on: ${{ matrix.os }} - strategy: - matrix: - os: [ubuntu-latest, macOS-latest, windows-latest] - steps: - - uses: actions/checkout@v4 - - uses: julia-actions/setup-julia@v1 - with: - version: nightly - arch: x64 - - uses: julia-actions/julia-buildpkg@latest - - uses: julia-actions/julia-runtest@latest - - uses: julia-actions/julia-uploadcodecov@latest - env: - CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }} diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 20622f0f..f0a6f0c3 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -11,11 +11,8 @@ jobs: runs-on: ${{ matrix.os }} strategy: matrix: - julia-version: ["1.0", "1.6", "1.9"] + julia-version: ["1.6", "1.9", "~1.10.0-0"] os: [ubuntu-latest, macOS-latest, windows-latest] - exclude: - - os: macOS-latest - julia-version: "1.0" steps: - uses: actions/checkout@v4 - uses: julia-actions/setup-julia@v1 diff --git a/NEWS.md b/NEWS.md index 5627adaa..ee9fae46 100644 --- a/NEWS.md +++ b/NEWS.md @@ -5,11 +5,37 @@ 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.15.0] xx/xx/2023 + +### Added + +- `ProductManifold` type was migrated from Manifolds.jl. +- `Fiber`, `VectorSpaceFiber` and `TangentSpace` types. `TangentSpace` is a generalized version of `TangentSpaceAtPoint` from Manifolds.jl. +- `change_representer!`, `change_metric!` and `Weingarten!` methods added to `PowerManifold`. +- `×` now also works for retractions, inverse retractions, and vector transports to create their product versions +- `retract`, `inverse_retract`, and `vector_transport_to` (and `_dir`) now also accept arbirtrary retractions on the product manifold. These act the same as the n-fold product of a retraction. + +### Changed + +- `Requires.jl` is added as a dependency to facilitate loading some methods related to `ProductManifolds` on Julia 1.6 to 1.8. Later versions rely on package extensions. +- `Documenter.jl` was updated to 1.0. +- `PowerManifold` can now store its size either in a field or in a type, similarly to `DefaultManifold`. By default the size is stored in a field. + +### Removed + +- Julia 1.0 is no longer supported. From now on, the earliest supported Julia version is 1.6. + +## [0.14.12] 23/09/2023 + +### Changed + +- Introduce a thorough way to allocate tangent vectors for `rand` + ## [0.14.11] 25/08/2023 ### Added -- MAke the `Weingarten` map a decorator capable function. +- Make the `Weingarten` map a decorator capable function. ## [0.14.10] 17/08/2023 diff --git a/Project.toml b/Project.toml index 33e92c23..2e4d8900 100644 --- a/Project.toml +++ b/Project.toml @@ -1,24 +1,34 @@ name = "ManifoldsBase" uuid = "3362f125-f0bb-47a3-aa74-596ffd7ef2fb" authors = ["Seth Axen ", "Mateusz Baran ", "Ronny Bergmann ", "Antoine Levitt "] -version = "0.14.12" +version = "0.15.0" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +Requires = "ae029012-a4dd-5104-9daa-d747884805df" + +[weakdeps] +RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" + +[extensions] +ManifoldsBaseRecursiveArrayToolsExt = "RecursiveArrayTools" [compat] DoubleFloats = ">= 0.9.2" +RecursiveArrayTools = "2" +Requires = "1" julia = "1.0" [extras] DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" +RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Test", "DoubleFloats", "ForwardDiff", "OrdinaryDiffEq", "ReverseDiff", "StaticArrays"] +test = ["Test", "DoubleFloats", "ForwardDiff", "OrdinaryDiffEq", "ReverseDiff", "StaticArrays", "RecursiveArrayTools"] diff --git a/docs/Project.toml b/docs/Project.toml index 9629fb2c..92399771 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,9 +1,11 @@ [deps] CondaPkg = "992eb4ea-22a4-4c89-a5bb-47a3300528ab" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244" ManifoldsBase = "3362f125-f0bb-47a3-aa74-596ffd7ef2fb" [compat] CondaPkg = "0.2" -Documenter = "0.27" -ManifoldsBase = "0.14" +Documenter = "1" +DocumenterCitations = "1.2" +ManifoldsBase = "0.15" diff --git a/docs/make.jl b/docs/make.jl index bdfc8d99..7a3e2f93 100755 --- a/docs/make.jl +++ b/docs/make.jl @@ -30,9 +30,13 @@ if "--quarto" ∈ ARGS end using Documenter: DocMeta, HTML, MathJax3, deploydocs, makedocs +using DocumenterCitations using ManifoldsBase -makedocs( +# (e) ...finally! make docs +bib = CitationBibliography(joinpath(@__DIR__, "src", "references.bib"); style = :alpha) +makedocs(; + # for development, we disable prettyurls format = HTML(; mathengine = MathJax3(), prettyurls = get(ENV, "CI", nothing) == "true", @@ -53,8 +57,11 @@ makedocs( "Vector transports" => "vector_transports.md", ], "Manifolds" => "manifolds.md", + "Meta-Manifolds" => "metamanifolds.md", "Decorating/Extending a Manifold" => "decorator.md", "Bases for tangent spaces" => "bases.md", + "References" => "references.md", ], + plugins = [bib], ) deploydocs(repo = "github.com/JuliaManifolds/ManifoldsBase.jl.git", push_preview = true) diff --git a/docs/src/assets/citations.css b/docs/src/assets/citations.css new file mode 100644 index 00000000..5ea61188 --- /dev/null +++ b/docs/src/assets/citations.css @@ -0,0 +1,19 @@ +/* Taken from https://juliadocs.org/DocumenterCitations.jl/v1.2/styling/ */ + +.citation dl { + display: grid; + grid-template-columns: max-content auto; } +.citation dt { + grid-column-start: 1; } +.citation dd { + grid-column-start: 2; + margin-bottom: 0.75em; } +.citation ul { + padding: 0 0 2.25em 0; + margin: 0; + list-style: none;} +.citation ul li { + text-indent: -2.25em; + margin: 0.33em 0.5em 0.5em 2.25em;} +.citation ol li { + padding-left:0.75em;} diff --git a/docs/src/manifolds.md b/docs/src/manifolds.md index 2f867ba2..78f6bab6 100644 --- a/docs/src/manifolds.md +++ b/docs/src/manifolds.md @@ -3,29 +3,7 @@ While the interface `ManifoldsBase.jl` does not cover concrete manifolds, it provides a few helpers to build or create manifolds based on existing manifolds -## [(Abstract) power manifold](@id sec-power-manifold) - -A power manifold is constructed like higher dimensional vector spaces are formed from the real line, just that for every point ``p = (p_1,\ldots,p_n) ∈ \mathcal M^n`` on the power manifold ``\mathcal M^n`` the entries of ``p`` are points ``p_1,\ldots,p_n ∈ \mathcal M`` on some manifold ``\mathcal M``. Note that ``n`` can also be replaced by multiple values, such that ``p`` is not a vector but a matrix or a multi-index array of points. - -```@autodocs -Modules = [ManifoldsBase] -Pages = ["src/PowerManifold.jl"] -Order = [:macro, :type, :function] -``` - -## `ValidationManifold` - -[`ValidationManifold`](@ref) is a simple decorator using the [`AbstractDecoratorManifold`](@ref) that “decorates” a manifold with tests that all involved points and vectors are valid for the wrapped manifold. -For example involved input and output paratemers are checked before and after running a function, repectively. -This is done by calling [`is_point`](@ref) or [`is_vector`](@ref) whenever applicable. - -```@autodocs -Modules = [ManifoldsBase] -Pages = ["ValidationManifold.jl"] -Order = [:macro, :type, :function] -``` - -## `DefaultManifold` +## A default manifold [`DefaultManifold`](@ref ManifoldsBase.DefaultManifold) is a simplified version of [`Euclidean`](https://juliamanifolds.github.io/Manifolds.jl/latest/manifolds/euclidean.html) and demonstrates a basic interface implementation. It can be used to perform simple tests. @@ -71,3 +49,16 @@ Modules = [ManifoldsBase] Pages = ["metric.jl"] Order = [:type, :function] ``` + + +## A manifold for validation + +[`ValidationManifold`](@ref) is a simple decorator using the [`AbstractDecoratorManifold`](@ref) that “decorates” a manifold with tests that all involved points and vectors are valid for the wrapped manifold. +For example involved input and output paratemers are checked before and after running a function, repectively. +This is done by calling [`is_point`](@ref) or [`is_vector`](@ref) whenever applicable. + +```@autodocs +Modules = [ManifoldsBase] +Pages = ["ValidationManifold.jl"] +Order = [:macro, :type, :function] +``` diff --git a/docs/src/metamanifolds.md b/docs/src/metamanifolds.md new file mode 100644 index 00000000..ad80e834 --- /dev/null +++ b/docs/src/metamanifolds.md @@ -0,0 +1,37 @@ +# Meta Manifolds + +While the interface does not provide concrete manifolds itself, it does provide several manifolds that can be build based on a given [`AbstractManifold`](@ref) instance. + +## [(Abstract) power manifold](@id sec-power-manifold) + +A power manifold is constructed like higher dimensional vector spaces are formed from the real line, just that for every point ``p = (p_1,\ldots,p_n) ∈ \mathcal M^n`` on the power manifold ``\mathcal M^n`` the entries of ``p`` are points ``p_1,\ldots,p_n ∈ \mathcal M`` on some manifold ``\mathcal M``. Note that ``n`` can also be replaced by multiple values, such that ``p`` is not a vector but a matrix or a multi-index array of points. + +```@autodocs +Modules = [ManifoldsBase] +Pages = ["src/PowerManifold.jl"] +Order = [:macro, :type, :function] +``` + +## [Product Manifold](@id ProductManifold) + +```@autodocs +Modules = [ManifoldsBase] +Pages = ["src/ProductManifold.jl"] +Order = [:macro, :type, :function] +``` + +## Fiber + +```@autodocs +Modules = [ManifoldsBase] +Pages = ["Fiber.jl"] +Order = [:macro, :type, :function] +``` + +## Tangent Space + +```@autodocs +Modules = [ManifoldsBase] +Pages = ["TangentSpace.jl"] +Order = [:macro, :type, :function] +``` diff --git a/docs/src/references.bib b/docs/src/references.bib new file mode 100644 index 00000000..a19a9d81 --- /dev/null +++ b/docs/src/references.bib @@ -0,0 +1,87 @@ +# ManifoldBase.jl +# +# Literature used within the Documentation of ManifoldsBase.jl +# ======================================================== +# +# citekeys should be of the form AllAuthors:Year, unless there is really many authors. +# +# ---------------------------------------------------------------------------------------- +@book{AbsilMahonySepulchre:2008, + AUTHOR = {Absil, P.-A. and Mahony, R. and Sepulchre, R.}, + DOI = {10.1515/9781400830244}, + NOTE = {available online at [press.princeton.edu/chapters/absil/](http://press.princeton.edu/chapters/absil/)}, + PUBLISHER = {Princeton University Press}, + TITLE = {Optimization Algorithms on Matrix Manifolds}, + YEAR = {2008}, +} + +@article{EhlersPiraniSchild:1972, + DOI = {10.1007/s10714-012-1353-4}, + YEAR = {1972}, + PUBLISHER = {Springer Science and Business Media {LLC}}, + VOLUME = {44}, + NUMBER = {6}, + PAGES = {1587--1609}, + AUTHOR = {J\"{u}rgen Ehlers and Felix A. E. Pirani and Alfred Schild}, + TITLE = {Republication of: The geometry of free fall and light propagation}, + JOURNAL = {General Relativity and Gravitation} +} + +@article{LorenziPennec:2013, + DOI = {10.1007/s10851-013-0470-3}, + EPRINT = {00870489}, + EPRINTTYPE = {HAL}, + YEAR = {2013}, + VOLUME = {50}, + NUMBER = {1-2}, + PAGES = {5--17}, + AUTHOR = {Marco Lorenzi and Xavier Pennec}, + TITLE = {Efficient Parallel Transport of Deformations in Time Series of Images: From Schild's to Pole Ladder}, + JOURNAL = {Journal of Mathematical Imaging and Vision} +} +@book{Lee:2019, + AUTHOR = {John M. Lee}, + DOI = {10.1007/978-3-319-91755-9}, + ISBN = {978-3-319-91755-9}, + PUBLISHER = {Springer Cham}, + TITLE = {Introduction to Riemannian Manifolds}, + YEAR = {2019} +} +@inproceedings{MuralidharanFletcher:2012, + DOI = {10.1109/cvpr.2012.6247780}, + YEAR = 2012, + AUTHOR = {P. Muralidharan and P. T. Fletcher}, + TITLE = {Sasaki metrics for analysis of longitudinal data on manifolds}, + BOOKTITLE = {2012 {IEEE} Conference on Computer Vision and Pattern Recognition} +} +@article{Pennec:2018, + AUTHOR ={Xavier Pennec}, + EPRINT = {1805.11436}, + EPRINTTYPE = {arXiv}, + Journal = {arXiv Preprint}, + TITLE = {Parallel Transport with Pole Ladder: a Third Order Scheme in Affine Connection Spaces which is Exact in Affine Symmetric Spaces.}, + YEAR = {2018}, + URL = {https://arxiv.org/abs/1805.11436}, +} +@article{Sasaki:1958, + DOI = {10.2748/tmj/1178244668}, + YEAR = {1958}, + VOLUME = {10}, + NUMBER = {3}, + AUTHOR = {Shigeo Sasaki}, + TITLE = {On the differential geometry of tangent bundles of Riemannian manifolds}, + JOURNAL = {Tohoku Math. J.} +} +@article{SatoIwai:2013, + AUTHOR = {Hiroyuki Sato and Toshihiro Iwai}, + ARCHIVEPREFIX = {arXiv}, + DOI = {10.1080/02331934.2013.836650}, + EPRINT = {1302.0125}, + JOURNAL = {Optimization}, + NUMBER = {4}, + PAGES = {1011--1031}, + PUBLISHER = {Informa {UK} Limited}, + TITLE = {A new, globally convergent Riemannian conjugate gradient method}, + VOLUME = {64}, + YEAR = {2013} +} \ No newline at end of file diff --git a/docs/src/references.md b/docs/src/references.md new file mode 100644 index 00000000..1bfafad1 --- /dev/null +++ b/docs/src/references.md @@ -0,0 +1,4 @@ +# Literature + +```@bibliography +``` \ No newline at end of file diff --git a/docs/src/types.md b/docs/src/types.md index ab94a674..d99468d5 100644 --- a/docs/src/types.md +++ b/docs/src/types.md @@ -9,7 +9,12 @@ Throughout the documentation of `ManifoldsBase.jl` we might use the [Euclidean S AbstractManifold ``` -which should store information about the manifold, for example parameters inherent to the manifold. +which should store information about the manifold, for example parameters inherent to the manifold. The parameters are stored in two possible ways, as a type parameter to dispatch on or as a field. For these the following internal functions exist + +```@docs +ManifoldsBase.wrap_type_parameter +ManifoldsBase.TypeParameter +``` ## Points on a manifold diff --git a/docs/src/vector_transports.md b/docs/src/vector_transports.md index 76e681c6..14936d01 100644 --- a/docs/src/vector_transports.md +++ b/docs/src/vector_transports.md @@ -7,7 +7,7 @@ A _vector transport_ is a way to transport a vector between two tangent spaces. Let ``p,q ∈ \mathcal M`` be given, ``c`` the curve along which we want to transport (cf. [parallel transport](@ref subsec-parallel-transport), for example a geodesic or curve given by a retraction. We can specify the geodesic or curve a retraction realises for example by a direction ``d``. -More precisely using [^AbsilMahonySepulchre2008], Def. 8.1.1, a vector transport +More precisely using [AbsilMahonySepulchre:2008](@cite), Def. 8.1.1, a vector transport ``T_{p,d}: T_p\mathcal M \to T_q\mathcal M``, ``p∈ \mathcal M``, ``Y∈ T_p\mathcal M`` is a smooth mapping associated to a retraction ``\operatorname{retr}_p(Y) = q`` such that diff --git a/ext/ManifoldsBaseRecursiveArrayToolsExt/ManifoldsBaseRecursiveArrayToolsExt.jl b/ext/ManifoldsBaseRecursiveArrayToolsExt/ManifoldsBaseRecursiveArrayToolsExt.jl new file mode 100644 index 00000000..3afeff9f --- /dev/null +++ b/ext/ManifoldsBaseRecursiveArrayToolsExt/ManifoldsBaseRecursiveArrayToolsExt.jl @@ -0,0 +1,69 @@ +module ManifoldsBaseRecursiveArrayToolsExt + +if isdefined(Base, :get_extension) + using ManifoldsBase + using RecursiveArrayTools + using ManifoldsBase: AbstractBasis, ProductBasisData, TangentSpaceType + using ManifoldsBase: number_of_components + using Random + + import ManifoldsBase: + allocate, + allocate_result, + default_inverse_retraction_method, + default_retraction_method, + default_vector_transport_method, + _get_dim_ranges, + get_vector, + get_vectors, + inverse_retract, + parallel_transport_direction, + parallel_transport_to, + project, + _retract, + riemann_tensor, + submanifold_component, + submanifold_components, + vector_transport_direction, + _vector_transport_direction, + _vector_transport_to, + vector_transport_to, + ziptuples + import Base: copyto!, getindex, setindex!, view +else + # imports need to be relative for Requires.jl-based workflows: + # https://github.com/JuliaArrays/ArrayInterface.jl/pull/387 + using ..ManifoldsBase + using ..RecursiveArrayTools + using ..ManifoldsBase: AbstractBasis, ProductBasisData, TangentSpaceType + using ..ManifoldsBase: number_of_components + using ..Random + + import ..ManifoldsBase: + allocate, + allocate_result, + default_inverse_retraction_method, + default_retraction_method, + default_vector_transport_method, + _get_dim_ranges, + get_vector, + get_vectors, + inverse_retract, + parallel_transport_direction, + parallel_transport_to, + project, + _retract, + riemann_tensor, + submanifold_component, + submanifold_components, + vector_transport_direction, + _vector_transport_direction, + _vector_transport_to, + vector_transport_to, + ziptuples + import ..Base: copyto!, getindex, setindex!, view +end + +include("ProductManifoldRecursiveArrayToolsExt.jl") + +end diff --git a/ext/ManifoldsBaseRecursiveArrayToolsExt/ProductManifoldRecursiveArrayToolsExt.jl b/ext/ManifoldsBaseRecursiveArrayToolsExt/ProductManifoldRecursiveArrayToolsExt.jl new file mode 100644 index 00000000..44367ffe --- /dev/null +++ b/ext/ManifoldsBaseRecursiveArrayToolsExt/ProductManifoldRecursiveArrayToolsExt.jl @@ -0,0 +1,375 @@ + +allocate(a::AbstractArray{<:ArrayPartition}) = map(allocate, a) +allocate(x::ArrayPartition) = ArrayPartition(map(allocate, x.x)...) +function allocate(x::ArrayPartition, T::Type) + return ArrayPartition(map(t -> allocate(t, T), submanifold_components(x))...) +end + +@inline function allocate_result(M::ProductManifold, f) + return ArrayPartition(map(N -> allocate_result(N, f), M.manifolds)) +end + +function copyto!(M::ProductManifold, q::ArrayPartition, p::ArrayPartition) + map(copyto!, M.manifolds, submanifold_components(q), submanifold_components(p)) + return q +end +function copyto!( + M::ProductManifold, + Y::ArrayPartition, + p::ArrayPartition, + X::ArrayPartition, +) + map( + copyto!, + M.manifolds, + submanifold_components(Y), + submanifold_components(p), + submanifold_components(X), + ) + return Y +end + + +function default_retraction_method(M::ProductManifold, ::Type{T}) where {T<:ArrayPartition} + return ProductRetraction( + map(default_retraction_method, M.manifolds, T.parameters[2].parameters)..., + ) +end + +function default_inverse_retraction_method( + M::ProductManifold, + ::Type{T}, +) where {T<:ArrayPartition} + return InverseProductRetraction( + map(default_inverse_retraction_method, M.manifolds, T.parameters[2].parameters)..., + ) +end + +function default_vector_transport_method( + M::ProductManifold, + ::Type{T}, +) where {T<:ArrayPartition} + return ProductVectorTransport( + map(default_vector_transport_method, M.manifolds, T.parameters[2].parameters)..., + ) +end + +function Base.exp(M::ProductManifold, p::ArrayPartition, X::ArrayPartition) + return ArrayPartition( + map( + exp, + M.manifolds, + submanifold_components(M, p), + submanifold_components(M, X), + )..., + ) +end +function Base.exp(M::ProductManifold, p::ArrayPartition, X::ArrayPartition, t::Number) + return ArrayPartition( + map( + (N, pc, Xc) -> exp(N, pc, Xc, t), + M.manifolds, + submanifold_components(M, p), + submanifold_components(M, X), + )..., + ) +end + +function get_vector( + M::ProductManifold, + p::ArrayPartition, + Xⁱ, + B::AbstractBasis{𝔽,TangentSpaceType}, +) where {𝔽} + dims = map(manifold_dimension, M.manifolds) + @assert length(Xⁱ) == sum(dims) + dim_ranges = _get_dim_ranges(dims) + tXⁱ = map(dr -> (@inbounds view(Xⁱ, dr)), dim_ranges) + ts = ziptuples(M.manifolds, submanifold_components(M, p), tXⁱ) + return ArrayPartition(map((@inline t -> get_vector(t..., B)), ts)) +end +function get_vector( + M::ProductManifold, + p::ArrayPartition, + Xⁱ, + B::CachedBasis{𝔽,<:AbstractBasis{𝔽},<:ProductBasisData}, +) where {𝔽} + dims = map(manifold_dimension, M.manifolds) + @assert length(Xⁱ) == sum(dims) + dim_ranges = _get_dim_ranges(dims) + tXⁱ = map(dr -> (@inbounds view(Xⁱ, dr)), dim_ranges) + ts = ziptuples(M.manifolds, submanifold_components(M, p), tXⁱ, B.data.parts) + return ArrayPartition(map((@inline t -> get_vector(t...)), ts)) +end + +function get_vectors( + M::ProductManifold, + p::ArrayPartition, + B::CachedBasis{𝔽,<:AbstractBasis{𝔽},<:ProductBasisData}, +) where {𝔽} + N = number_of_components(M) + xparts = submanifold_components(p) + BVs = map(t -> get_vectors(t...), ziptuples(M.manifolds, xparts, B.data.parts)) + zero_tvs = map(t -> zero_vector(t...), ziptuples(M.manifolds, xparts)) + vs = typeof(ArrayPartition(zero_tvs...))[] + for i in 1:N, k in 1:length(BVs[i]) + push!( + vs, + ArrayPartition(zero_tvs[1:(i - 1)]..., BVs[i][k], zero_tvs[(i + 1):end]...), + ) + end + return vs +end + +ManifoldsBase._get_vector_cache_broadcast(::ArrayPartition) = Val(false) + +""" + getindex(p, M::ProductManifold, i::Union{Integer,Colon,AbstractVector}) + p[M::ProductManifold, i] + +Access the element(s) at index `i` of a point `p` on a [`ProductManifold`](@ref) `M` by +linear indexing. +See also [Array Indexing](https://docs.julialang.org/en/v1/manual/arrays/#man-array-indexing-1) in Julia. +""" +Base.@propagate_inbounds function Base.getindex( + p::ArrayPartition, + M::ProductManifold, + i::Union{Integer,Colon,AbstractVector,Val}, +) + return get_component(M, p, i) +end + +function inverse_retract( + M::ProductManifold, + p::ArrayPartition, + q::ArrayPartition, + method::InverseProductRetraction, +) + return ArrayPartition( + map( + inverse_retract, + M.manifolds, + submanifold_components(M, p), + submanifold_components(M, q), + method.inverse_retractions, + ), + ) +end + +function Base.log(M::ProductManifold, p::ArrayPartition, q::ArrayPartition) + return ArrayPartition( + map( + log, + M.manifolds, + submanifold_components(M, p), + submanifold_components(M, q), + )..., + ) +end + +function parallel_transport_direction( + M::ProductManifold, + p::ArrayPartition, + X::ArrayPartition, + d::ArrayPartition, +) + return ArrayPartition( + map( + parallel_transport_direction, + M.manifolds, + submanifold_components(M, p), + submanifold_components(M, X), + submanifold_components(M, d), + ), + ) +end + +function parallel_transport_to( + M::ProductManifold, + p::ArrayPartition, + X::ArrayPartition, + q::ArrayPartition, +) + return ArrayPartition( + map( + parallel_transport_to, + M.manifolds, + submanifold_components(M, p), + submanifold_components(M, X), + submanifold_components(M, q), + ), + ) +end + +function project(M::ProductManifold, p::ArrayPartition) + return ArrayPartition(map(project, M.manifolds, submanifold_components(M, p))...) +end +function project(M::ProductManifold, p::ArrayPartition, X::ArrayPartition) + return ArrayPartition( + map( + project, + M.manifolds, + submanifold_components(M, p), + submanifold_components(M, X), + )..., + ) +end + +@doc raw""" + rand(M::ProductManifold; parts_kwargs = map(_ -> (;), M.manifolds)) + +Return a random point on [`ProductManifold`](@ref) `M`. `parts_kwargs` is +a tuple of keyword arguments for `rand` on each manifold in `M.manifolds`. +""" +function Random.rand( + M::ProductManifold; + vector_at = nothing, + parts_kwargs = map(_ -> (;), M.manifolds), +) + if vector_at === nothing + return ArrayPartition( + map((N, kwargs) -> rand(N; kwargs...), M.manifolds, parts_kwargs)..., + ) + else + return ArrayPartition( + map( + (N, p, kwargs) -> rand(N; vector_at = p, kwargs...), + M.manifolds, + submanifold_components(M, vector_at), + parts_kwargs, + )..., + ) + end +end +function Random.rand( + rng::AbstractRNG, + M::ProductManifold; + vector_at = nothing, + parts_kwargs = map(_ -> (;), M.manifolds), +) + if vector_at === nothing + return ArrayPartition( + map((N, kwargs) -> rand(rng, N; kwargs...), M.manifolds, parts_kwargs)..., + ) + else + return ArrayPartition( + map( + (N, p, kwargs) -> rand(rng, N; vector_at = p, kwargs...), + M.manifolds, + submanifold_components(M, vector_at), + parts_kwargs, + )..., + ) + end +end + +function _retract( + M::ProductManifold, + p::ArrayPartition, + X::ArrayPartition, + t::Number, + method::ProductRetraction, +) + return ArrayPartition( + map( + (N, pc, Xc, rm) -> retract(N, pc, Xc, t, rm), + M.manifolds, + submanifold_components(M, p), + submanifold_components(M, X), + method.retractions, + ), + ) +end + +function riemann_tensor( + M::ProductManifold, + p::ArrayPartition, + X::ArrayPartition, + Y::ArrayPartition, + Z::ArrayPartition, +) + return ArrayPartition( + map( + riemann_tensor, + M.manifolds, + submanifold_components(M, p), + submanifold_components(M, X), + submanifold_components(M, Y), + submanifold_components(M, Z), + ), + ) +end + +""" + setindex!(q, p, M::ProductManifold, i::Union{Integer,Colon,AbstractVector}) + q[M::ProductManifold,i...] = p + +set the element `[i...]` of a point `q` on a [`ProductManifold`](@ref) by linear indexing to `q`. +See also [Array Indexing](https://docs.julialang.org/en/v1/manual/arrays/#man-array-indexing-1) in Julia. +""" +Base.@propagate_inbounds function Base.setindex!( + q::ArrayPartition, + p, + M::ProductManifold, + i::Union{Integer,Colon,AbstractVector,Val}, +) + return set_component!(M, q, p, i) +end + +@inline submanifold_component(p::ArrayPartition, ::Val{I}) where {I} = p.x[I] +@inline submanifold_components(p::ArrayPartition) = p.x + +function vector_transport_direction( + M::ProductManifold, + p::ArrayPartition, + X::ArrayPartition, + d::ArrayPartition, + m::ProductVectorTransport, +) + return ArrayPartition( + map( + vector_transport_direction, + M.manifolds, + submanifold_components(M, p), + submanifold_components(M, X), + submanifold_components(M, d), + m.methods, + ), + ) +end + +function vector_transport_to( + M::ProductManifold, + p::ArrayPartition, + X::ArrayPartition, + q::ArrayPartition, + m::ProductVectorTransport, +) + return ArrayPartition( + map( + vector_transport_to, + M.manifolds, + submanifold_components(M, p), + submanifold_components(M, X), + submanifold_components(M, q), + m.methods, + ), + ) +end +function vector_transport_to( + M::ProductManifold, + p::ArrayPartition, + X::ArrayPartition, + q::ArrayPartition, + m::ParallelTransport, +) + return ArrayPartition( + map( + (iM, ip, iX, id) -> vector_transport_to(iM, ip, iX, id, m), + M.manifolds, + submanifold_components(M, p), + submanifold_components(M, X), + submanifold_components(M, q), + ), + ) +end diff --git a/src/DefaultManifold.jl b/src/DefaultManifold.jl index 2ff34103..c10957c3 100644 --- a/src/DefaultManifold.jl +++ b/src/DefaultManifold.jl @@ -29,6 +29,14 @@ function DefaultManifold(n::Vararg{Int}; field = ℝ, parameter::Symbol = :field return DefaultManifold{field,typeof(size)}(size) end +function allocation_promotion_function( + ::DefaultManifold{ℂ}, + ::Union{typeof(get_vector),typeof(get_coordinates)}, + ::Tuple, +) + return complex +end + change_representer!(M::DefaultManifold, Y, ::EuclideanMetric, p, X) = copyto!(M, Y, p, X) change_metric!(M::DefaultManifold, Y, ::EuclideanMetric, p, X) = copyto!(M, Y, p, X) diff --git a/src/Fiber.jl b/src/Fiber.jl new file mode 100644 index 00000000..5b4feab6 --- /dev/null +++ b/src/Fiber.jl @@ -0,0 +1,51 @@ + +""" + abstract type FiberType end + +An abstract type for fiber types that can be used within [`Fiber`](@ref). +""" +abstract type FiberType end + +@doc raw""" + Fiber{𝔽,TFiber<:FiberType,TM<:AbstractManifold{𝔽},TX} <: AbstractManifold{𝔽} + +A fiber of a fiber bundle at a point `p` on the manifold. + +This fiber itself is also a `manifold`. For vector fibers it's by default flat and hence +isometric to the [`Euclidean`](https://juliamanifolds.github.io/Manifolds.jl/latest/manifolds/euclidean.html) manifold. + +# Fields + +* `manifold` – base space of the fiber bundle +* `point` – a point ``p`` from the base space; the fiber corresponds to the preimage + by bundle projection ``\pi^{-1}(\{p\})``. + + +# Constructor + + Fiber(M::AbstractManifold, p, fiber_type::FiberType) + +A fiber of type `fiber_type` at point `p` from the manifold `manifold`. +""" +struct Fiber{𝔽,TFiber<:FiberType,TM<:AbstractManifold{𝔽},TX} <: AbstractManifold{𝔽} + manifold::TM + point::TX + fiber_type::TFiber +end + +base_manifold(B::Fiber) = B.manifold + +function Base.show(io::IO, ::MIME"text/plain", vs::Fiber) + summary(io, vs) + println(io, "\nFiber:") + pre = " " + sf = sprint(show, "text/plain", vs.fiber_type; context = io, sizehint = 0) + sf = replace(sf, '\n' => "\n$(pre)") + sm = sprint(show, "text/plain", vs.manifold; context = io, sizehint = 0) + sm = replace(sm, '\n' => "\n$(pre)") + println(io, pre, sf, sm) + println(io, "Base point:") + sp = sprint(show, "text/plain", vs.point; context = io, sizehint = 0) + sp = replace(sp, '\n' => "\n$(pre)") + return print(io, pre, sp) +end diff --git a/src/ManifoldsBase.jl b/src/ManifoldsBase.jl index 9296f1fc..36707dc1 100644 --- a/src/ManifoldsBase.jl +++ b/src/ManifoldsBase.jl @@ -17,16 +17,17 @@ import Base: -, *, == -import LinearAlgebra: dot, norm, det, cross, I, UniformScaling, Diagonal - +import LinearAlgebra: ×, dot, norm, det, cross, I, UniformScaling, Diagonal import Random: rand, rand! -import Markdown: @doc_str using LinearAlgebra +using Markdown: @doc_str using Random +using Requires include("maintypes.jl") include("numbers.jl") +include("Fiber.jl") include("bases.jl") include("retractions.jl") include("exp_log_geo.jl") @@ -797,18 +798,6 @@ function mid_point!(M::AbstractManifold, q, p1, p2) return exp!(M, q, p1, X / 2) end -@static if VERSION <= v"1.1" - function mid_point!( - M::AbstractManifold, - q::AbstractArray{T1,0}, - p1::AbstractArray{T2,0}, - p2::AbstractArray{T3,0}, - ) where {T1,T2,T3} - X = log(M, p1, p2) - return exp!(M, q, p1, fill(X / 2)) - end -end - """ norm(M::AbstractManifold, p, X) @@ -978,12 +967,38 @@ include("vector_spaces.jl") include("point_vector_fallbacks.jl") include("nested_trait.jl") include("decorator_trait.jl") + +include("VectorFiber.jl") +include("TangentSpace.jl") include("ValidationManifold.jl") include("EmbeddedManifold.jl") include("DefaultManifold.jl") +include("ProductManifold.jl") include("PowerManifold.jl") +# +# +# Requires +# ----- +function __init__() + @static if !isdefined(Base, :get_extension) + @require RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" begin + include( + "../ext/ManifoldsBaseRecursiveArrayToolsExt/ManifoldsBaseRecursiveArrayToolsExt.jl", + ) + end + end +end +# +# +# Export +# ------ +# +# (a) Manifolds and general types export AbstractManifold, AbstractManifoldPoint, TVector, CoTVector, TFVector, CoTFVector +export VectorSpaceFiber +export TangentSpace, TangentSpaceType +export CotangentSpace, CotangentSpaceType export AbstractDecoratorManifold export AbstractTrait, IsEmbeddedManifold, IsEmbeddedSubmanifold, IsIsometricEmbeddedManifold export IsExplicitDecorator @@ -992,9 +1007,9 @@ export EmbeddedManifold export AbstractPowerManifold, PowerManifold export AbstractPowerRepresentation, NestedPowerRepresentation, NestedReplacingPowerRepresentation +export ProductManifold -export OutOfInjectivityRadiusError, ManifoldDomainError - +# (b) Retraction Types export AbstractRetractionMethod, ApproximateInverseRetraction, CayleyRetraction, @@ -1005,15 +1020,19 @@ export AbstractRetractionMethod, QRRetraction, PadeRetraction, PolarRetraction, + ProductRetraction, ProjectionRetraction, RetractionWithKeywords, + SasakiRetraction, ShootingInverseRetraction, SoftmaxRetraction +# (c) Inverse Retraction Types export AbstractInverseRetractionMethod, ApproximateInverseRetraction, CayleyInverseRetraction, EmbeddedInverseRetraction, + InverseProductRetraction, LogarithmicInverseRetraction, NLSolveInverseRetraction, QRInverseRetraction, @@ -1023,10 +1042,12 @@ export AbstractInverseRetractionMethod, InverseRetractionWithKeywords, SoftmaxInverseRetraction +# (d) Vector Transport Types export AbstractVectorTransportMethod, DifferentiatedRetractionVectorTransport, ParallelTransport, PoleLadderTransport, + ProductVectorTransport, ProjectionTransport, ScaledVectorTransport, SchildsLadderTransport, @@ -1034,6 +1055,7 @@ export AbstractVectorTransportMethod, VectorTransportTo, VectorTransportWithKeywords +# (e) Basis Types export CachedBasis, DefaultBasis, DefaultOrthogonalBasis, @@ -1044,10 +1066,16 @@ export CachedBasis, ProjectedOrthonormalBasis, VeeOrthogonalBasis +# (f) Error Messages +export OutOfInjectivityRadiusError, ManifoldDomainError export ApproximatelyError export CompositeManifoldError, ComponentManifoldError, ManifoldDomainError -export allocate, +# (g) Functions on Manifolds +export ×, + ℝ, + ℂ, + allocate, angle, base_manifold, change_basis, @@ -1120,6 +1148,7 @@ export allocate, retract!, riemann_tensor, riemann_tensor!, + vector_space_dimension, vector_transport_along, vector_transport_along!, vector_transport_direction, @@ -1132,5 +1161,4 @@ export allocate, Weingarten!, zero_vector, zero_vector! - end # module diff --git a/src/PowerManifold.jl b/src/PowerManifold.jl index 61e45e0a..ec65a984 100644 --- a/src/PowerManifold.jl +++ b/src/PowerManifold.jl @@ -43,73 +43,94 @@ abstract type AbstractPowerManifold{ } <: AbstractManifold{𝔽} end @doc raw""" - PowerManifold{𝔽,TM<:AbstractManifold,TSize<:Tuple,TPR<:AbstractPowerRepresentation} <: AbstractPowerManifold{𝔽,TM} + PowerManifold{𝔽,TM<:AbstractManifold,TSize,TPR<:AbstractPowerRepresentation} <: AbstractPowerManifold{𝔽,TM} -The power manifold $\mathcal M^{n_1× n_2 × … × n_d}$ with power geometry - `TSize` statically defines the number of elements along each axis. +The power manifold ``\mathcal M^{n_1× n_2 × … × n_d}`` with power geometry. + `TSize` defines the number of elements along each axis, either statically using + `TypeParameter` or storing it in a field. For example, a manifold-valued time series would be represented by a power manifold with -$d$ equal to 1 and $n_1$ equal to the number of samples. A manifold-valued image +``d`` equal to 1 and ``n_1`` equal to the number of samples. A manifold-valued image (for example in diffusion tensor imaging) would be represented by a two-axis power -manifold ($d=2$) with $n_1$ and $n_2$ equal to width and height of the image. +manifold (``d=2``) with ``n_1`` and ``n_2`` equal to width and height of the image. While the size of the manifold is static, points on the power manifold would not be represented by statically-sized arrays. # Constructor - PowerManifold(M::PowerManifold, N_1, N_2, ..., N_d) - PowerManifold(M::AbstractManifold, NestedPowerRepresentation(), N_1, N_2, ..., N_d) + PowerManifold(M::PowerManifold, N_1, N_2, ..., N_d; parameter::Symbol=:field) + PowerManifold(M::AbstractManifold, NestedPowerRepresentation(), N_1, N_2, ..., N_d; parameter::Symbol=:field) M^(N_1, N_2, ..., N_d) -Generate the power manifold $M^{N_1 × N_2 × … × N_d}$. -By default, a [`PowerManifold`](@ref} is expanded further, i.e. for `M=PowerManifold(N,3)` -`PowerManifold(M,2)` is equivalend to `PowerManifold(N,3,2)`. Points are then 3×2 matrices +Generate the power manifold ``M^{N_1 × N_2 × … × N_d}``. +By default, a [`PowerManifold`](@ref) is expanded further, i.e. for `M=PowerManifold(N, 3)` +`PowerManifold(M, 2)` is equivalent to `PowerManifold(N, 3, 2)`. Points are then 3×2 matrices of points on `N`. Providing a [`NestedPowerRepresentation`](@ref) as the second argument to the constructor -can be used to nest manifold, i.e. `PowerManifold(M,NestedPowerRepresentation(),2)` +can be used to nest manifold, i.e. `PowerManifold(M, NestedPowerRepresentation(), 2)` represents vectors of length 2 whose elements are vectors of length 3 of points on N in a nested array representation. Since there is no default [`AbstractPowerRepresentation`](@ref) within this interface, the `^` operator is only available for `PowerManifold`s and concatenates dimensions. + +`parameter`: whether a type parameter should be used to store `n`. By default size +is stored in a field. Value can either be `:field` or `:type`. """ struct PowerManifold{𝔽,TM<:AbstractManifold{𝔽},TSize,TPR<:AbstractPowerRepresentation} <: AbstractPowerManifold{𝔽,TM,TPR} manifold::TM + size::TSize +end + +""" + _parameter_symbol(M::PowerManifold) + +Return `:field` if size of [`PowerManifold`](@ref) `M` is stored in a field and `:type` +if in a `TypeParameter`. +""" +_parameter_symbol(::PowerManifold) = :field +function _parameter_symbol( + ::PowerManifold{𝔽,<:AbstractManifold{𝔽},<:TypeParameter}, +) where {𝔽} + return :type end function PowerManifold( M::AbstractManifold{𝔽}, ::TPR, - size::Integer..., + size::Integer...; + parameter::Symbol = :field, ) where {𝔽,TPR<:AbstractPowerRepresentation} - return PowerManifold{𝔽,typeof(M),Tuple{size...},TPR}(M) + size_w = wrap_type_parameter(parameter, size) + return PowerManifold{𝔽,typeof(M),typeof(size_w),TPR}(M, size_w) end function PowerManifold( M::PowerManifold{𝔽,TM,TSize,TPR}, - size::Integer..., + size::Integer...; + parameter::Symbol = _parameter_symbol(M), ) where {𝔽,TM<:AbstractManifold{𝔽},TSize,TPR<:AbstractPowerRepresentation} - return PowerManifold{𝔽,TM,Tuple{TSize.parameters...,size...},TPR}(M.manifold) + size_w = wrap_type_parameter(parameter, (get_parameter(M.size)..., size...)) + return PowerManifold{𝔽,TM,typeof(size_w),TPR}(M.manifold, size_w) end function PowerManifold( - M::PowerManifold{𝔽,TM,TSize}, + M::PowerManifold{𝔽,TM}, ::TPR, - size::Integer..., -) where {𝔽,TM<:AbstractManifold{𝔽},TSize,TPR<:AbstractPowerRepresentation} - return PowerManifold{𝔽,TM,Tuple{TSize.parameters...,size...},TPR}(M.manifold) + size::Integer...; + parameter::Symbol = _parameter_symbol(M), +) where {𝔽,TM<:AbstractManifold{𝔽},TPR<:AbstractPowerRepresentation} + size_w = wrap_type_parameter(parameter, (get_parameter(M.size)..., size...)) + return PowerManifold{𝔽,TM,typeof(size_w),TPR}(M.manifold, size_w) end function PowerManifold( - M::PowerManifold{𝔽,TM,TSize}, + M::PowerManifold{𝔽}, ::TPR, - size::Integer..., -) where { - 𝔽, - TM<:AbstractManifold{𝔽}, - TSize, - TPR<:Union{NestedPowerRepresentation,NestedReplacingPowerRepresentation}, -} - return PowerManifold{𝔽,PowerManifold{𝔽,TM,TSize},Tuple{size...},TPR}(M) + size::Integer...; + parameter::Symbol = _parameter_symbol(M), +) where {𝔽,TPR<:Union{NestedPowerRepresentation,NestedReplacingPowerRepresentation}} + size_w = wrap_type_parameter(parameter, size) + return PowerManifold{𝔽,typeof(M),typeof(size_w),TPR}(M, size_w) end """ @@ -225,6 +246,47 @@ function allocation_promotion_function(M::AbstractPowerManifold, f, args::Tuple) return allocation_promotion_function(M.manifold, f, args) end +""" + change_representer(M::AbstractPowerManifold, ::AbstractMetric, p, X) + +Since the metric on a power manifold decouples, the change of a representer can be done elementwise +""" +change_representer(::AbstractPowerManifold, ::AbstractMetric, ::Any, ::Any) + +function change_representer!(M::AbstractPowerManifold, Y, G::AbstractMetric, p, X) + rep_size = representation_size(M.manifold) + for i in get_iterator(M) + change_representer!( + M.manifold, + _write(M, rep_size, Y, i), + G, + _read(M, rep_size, p, i), + _read(M, rep_size, X, i), + ) + end + return Y +end + +""" + change_metric(M::AbstractPowerManifold, ::AbstractMetric, p, X) + +Since the metric on a power manifold decouples, the change of metric can be done elementwise. +""" +change_metric(M::AbstractPowerManifold, ::AbstractMetric, ::Any, ::Any) + +function change_metric!(M::AbstractPowerManifold, Y, G::AbstractMetric, p, X) + rep_size = representation_size(M.manifold) + for i in get_iterator(M) + change_metric!( + M.manifold, + _write(M, rep_size, Y, i), + G, + _read(M, rep_size, p, i), + _read(M, rep_size, X, i), + ) + end + return Y +end """ check_point(M::AbstractPowerManifold, p; kwargs...) @@ -253,7 +315,8 @@ end check_power_size(M, p) check_power_size(M, p, X) -Check whether p hase the right size to represent points on M generically, i.e. just cheking the overall sizes, not the individual ones per manifold +Check whether `p`` has the right size to represent points on `M`` generically, i.e. just +checking the overall sizes, not the individual ones per manifold. """ function check_power_size(M::AbstractPowerManifold, p) d = prod(representation_size(M.manifold)) * prod(power_dimensions(M)) @@ -555,13 +618,24 @@ function get_coordinates!( return c end -get_iterator(::PowerManifold{𝔽,<:AbstractManifold{𝔽},Tuple{N}}) where {𝔽,N} = Base.OneTo(N) +function get_iterator( + ::PowerManifold{𝔽,<:AbstractManifold{𝔽},TypeParameter{Tuple{N}}}, +) where {𝔽,N} + return Base.OneTo(N) +end +function get_iterator(M::PowerManifold{𝔽,<:AbstractManifold{𝔽},Tuple{Int}}) where {𝔽} + return Base.OneTo(M.size[1]) +end @generated function get_iterator( - ::PowerManifold{𝔽,<:AbstractManifold{𝔽},SizeTuple}, + ::PowerManifold{𝔽,<:AbstractManifold{𝔽},TypeParameter{SizeTuple}}, ) where {𝔽,SizeTuple} size_tuple = size_to_tuple(SizeTuple) return Base.product(map(Base.OneTo, size_tuple)...) end +function get_iterator(M::PowerManifold{𝔽,<:AbstractManifold{𝔽},NTuple{N,Int}}) where {𝔽,N} + size_tuple = M.size + return Base.product(map(Base.OneTo, size_tuple)...) +end function get_vector( M::AbstractPowerManifold, @@ -880,16 +954,17 @@ end manifold_dimension(M::PowerManifold) Returns the manifold-dimension of an [`PowerManifold`](@ref) `M` -$=\mathcal N = (\mathcal M)^{n_1,…,n_d}$, i.e. with $n=(n_1,…,n_d)$ the array -size of the power manifold and $d_{\mathcal M}$ the dimension of the base manifold -$\mathcal M$, the manifold is of dimension +``=\mathcal N = (\mathcal M)^{n_1,…,n_d}``, i.e. with ``n=(n_1,…,n_d)`` the array +size of the power manifold and ``d_{\mathcal M}`` the dimension of the base manifold +``\mathcal M``, the manifold is of dimension ````math \dim(\mathcal N) = \dim(\mathcal M)\prod_{i=1}^d n_i = n_1n_2\cdot…\cdot n_d \dim(\mathcal M). ```` """ -function manifold_dimension(M::PowerManifold{𝔽,<:AbstractManifold,TSize}) where {𝔽,TSize} - return manifold_dimension(M.manifold) * prod(size_to_tuple(TSize)) +function manifold_dimension(M::PowerManifold) + size = get_parameter(M.size) + return manifold_dimension(M.manifold) * prod(size) end function mid_point!(M::AbstractPowerManifold, q, p1, p2) @@ -1010,8 +1085,8 @@ end return the power of `M`, """ -function power_dimensions(::PowerManifold{𝔽,<:AbstractManifold,TSize}) where {𝔽,TSize} - return size_to_tuple(TSize) +function power_dimensions(M::PowerManifold) + return get_parameter(M.size) end @doc raw""" @@ -1280,10 +1355,21 @@ Base.@propagate_inbounds function Base.setindex!( return set_component!(M, q, p, I...) end -function Base.show(io::IO, M::PowerManifold{𝔽,TM,TSize,TPR}) where {𝔽,TM,TSize,TPR} +function Base.show( + io::IO, + M::PowerManifold{𝔽,TM,TSize,TPR}, +) where {𝔽,TM<:AbstractManifold{𝔽},TSize,TPR<:AbstractPowerRepresentation} + size = get_parameter(M.size) + return print(io, "PowerManifold($(M.manifold), $(TPR()), $(join(size, ", ")))") +end +function Base.show( + io::IO, + M::PowerManifold{𝔽,TM,TypeParameter{TSize},TPR}, +) where {𝔽,TM<:AbstractManifold{𝔽},TSize,TPR<:AbstractPowerRepresentation} + size = get_parameter(M.size) return print( io, - "PowerManifold($(M.manifold), $(TPR()), $(join(TSize.parameters, ", ")))", + "PowerManifold($(M.manifold), $(TPR()), $(join(size, ", ")); parameter=:type)", ) end @@ -1456,6 +1542,29 @@ function Base.view( return view(p[I...], rep_size_to_colons(rep_size)...) end +@doc raw""" + Y = Weingarten(M::AbstractPowerManifold, p, X, V) + Weingarten!(M::AbstractPowerManifold, Y, p, X, V) + +Since the metric decouples, also the computation of the Weingarten map +``\mathcal W_p`` can be computed elementwise on the single elements of the [`PowerManifold`](@ref) `M`. +""" +Weingarten(::AbstractPowerManifold, p, X, V) + +function Weingarten!(M::AbstractPowerManifold, Y, p, X, V) + rep_size = representation_size(M.manifold) + for i in get_iterator(M) + Weingarten!( + M.manifold, + _write(M, rep_size, Y, i), + _read(M, rep_size, p, i), + _read(M, rep_size, X, i), + _read(M, rep_size, V, i), + ) + end + return Y +end + @inline function _write(M::AbstractPowerManifold, rep_size::Tuple, x::AbstractArray, i::Int) return _write(M, rep_size, x, (i,)) end diff --git a/src/ProductManifold.jl b/src/ProductManifold.jl new file mode 100644 index 00000000..a11983fe --- /dev/null +++ b/src/ProductManifold.jl @@ -0,0 +1,1212 @@ +@doc raw""" + ProductManifold{𝔽,TM<:Tuple} <: AbstractManifold{𝔽} + +Product manifold $M_1 × M_2 × … × M_n$ with product geometry. + +# Constructor + + ProductManifold(M_1, M_2, ..., M_n) + +generates the product manifold $M_1 × M_2 × … × M_n$. +Alternatively, the same manifold can be contructed using the `×` operator: +`M_1 × M_2 × M_3`. +""" +struct ProductManifold{𝔽,TM<:Tuple} <: AbstractDecoratorManifold{𝔽} + manifolds::TM +end + +function ProductManifold(manifolds::AbstractManifold...) + 𝔽 = ManifoldsBase._unify_number_systems((number_system.(manifolds))...) + return ProductManifold{𝔽,typeof(manifolds)}(manifolds) +end + +""" + getindex(M::ProductManifold, i) + M[i] + +access the `i`th manifold component from the [`ProductManifold`](@ref) `M`. +""" +@inline Base.getindex(M::ProductManifold, i::Integer) = M.manifolds[i] + +ProductManifold() = throw(MethodError("No method matching ProductManifold().")) + +const PRODUCT_BASIS_LIST = [ + VeeOrthogonalBasis, + DefaultBasis, + DefaultBasis{<:Any,TangentSpaceType}, + DefaultOrthogonalBasis, + DefaultOrthogonalBasis{<:Any,TangentSpaceType}, + DefaultOrthonormalBasis, + DefaultOrthonormalBasis{<:Any,TangentSpaceType}, + ProjectedOrthonormalBasis{:gram_schmidt,ℝ}, + ProjectedOrthonormalBasis{:svd,ℝ}, +] + +""" + ProductBasisData + +A typed tuple to store tuples of data of stored/precomputed bases for a [`ProductManifold`](@ref). +""" +struct ProductBasisData{T<:Tuple} + parts::T +end + +const PRODUCT_BASIS_LIST_CACHED = [CachedBasis] + +""" + ProductMetric <: AbstractMetric + +A type to represent the product of metrics for a [`ProductManifold`](@ref). +""" +struct ProductMetric <: AbstractMetric end + +""" + ProductRetraction(retractions::AbstractRetractionMethod...) + +Product retraction of `retractions`. Works on [`ProductManifold`](@ref). +""" +struct ProductRetraction{TR<:Tuple} <: AbstractRetractionMethod + retractions::TR +end + +function ProductRetraction(retractions::AbstractRetractionMethod...) + return ProductRetraction{typeof(retractions)}(retractions) +end + +""" + InverseProductRetraction(retractions::AbstractInverseRetractionMethod...) + +Product inverse retraction of `inverse retractions`. Works on [`ProductManifold`](@ref). +""" +struct InverseProductRetraction{TR<:Tuple} <: AbstractInverseRetractionMethod + inverse_retractions::TR +end + +function InverseProductRetraction(inverse_retractions::AbstractInverseRetractionMethod...) + return InverseProductRetraction{typeof(inverse_retractions)}(inverse_retractions) +end + +function allocation_promotion_function(M::ProductManifold, f, args::Tuple) + apfs = map(MM -> allocation_promotion_function(MM, f, args), M.manifolds) + return reduce(combine_allocation_promotion_functions, apfs) +end + +""" + ProductVectorTransport(methods::AbstractVectorTransportMethod...) + +Product vector transport type of `methods`. Works on [`ProductManifold`](@ref). +""" +struct ProductVectorTransport{TR<:Tuple} <: AbstractVectorTransportMethod + methods::TR +end + +function ProductVectorTransport(methods::AbstractVectorTransportMethod...) + return ProductVectorTransport{typeof(methods)}(methods) +end + +""" + change_representer(M::ProductManifold, ::AbstractMetric, p, X) + +Since the metric on a product manifold decouples, the change of a representer can be done elementwise +""" +change_representer(::ProductManifold, ::AbstractMetric, ::Any, ::Any) + +function change_representer!(M::ProductManifold, Y, G::AbstractMetric, p, X) + map( + (m, y, P, x) -> change_representer!(m, y, G, P, x), + M.manifolds, + submanifold_components(M, Y), + submanifold_components(M, p), + submanifold_components(M, X), + ) + return Y +end + +""" + change_metric(M::ProductManifold, ::AbstractMetric, p, X) + +Since the metric on a product manifold decouples, the change of metric can be done elementwise. +""" +change_metric(::ProductManifold, ::AbstractMetric, ::Any, ::Any) + +function change_metric!(M::ProductManifold, Y, G::AbstractMetric, p, X) + map( + (m, y, P, x) -> change_metric!(m, y, G, P, x), + M.manifolds, + submanifold_components(M, Y), + submanifold_components(M, p), + submanifold_components(M, X), + ) + return Y +end + +""" + check_point(M::ProductManifold, p; kwargs...) + +Check whether `p` is a valid point on the [`ProductManifold`](@ref) `M`. +If `p` is not a point on `M` a [`CompositeManifoldError`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/functions.html#ManifoldsBase.CompositeManifoldError).consisting of all error messages of the +components, for which the tests fail is returned. + +The tolerance for the last test can be set using the `kwargs...`. +""" +function check_point(M::ProductManifold, p; kwargs...) + try + submanifold_components(M, p) + catch e + return DomainError("Point $p does not support submanifold_components") + end + ts = ziptuples(Tuple(1:length(M.manifolds)), M.manifolds, submanifold_components(M, p)) + e = [(t[1], check_point(t[2:end]...; kwargs...)) for t in ts] + errors = filter((x) -> !(x[2] === nothing), e) + cerr = [ComponentManifoldError(er...) for er in errors] + (length(errors) > 1) && return CompositeManifoldError(cerr) + (length(errors) == 1) && return cerr[1] + return nothing +end + +""" + check_size(M::ProductManifold, p; kwargs...) + +Check whether `p` is of valid size on the [`ProductManifold`](@ref) `M`. +If `p` has components of wrong size a [`CompositeManifoldError`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/functions.html#ManifoldsBase.CompositeManifoldError).consisting of all error messages of the +components, for which the tests fail is returned. + +The tolerance for the last test can be set using the `kwargs...`. +""" +function check_size(M::ProductManifold, p) + try + submanifold_components(M, p) + catch e + return DomainError("Point $p does not support submanifold_components") + end + ts = ziptuples(Tuple(1:length(M.manifolds)), M.manifolds, submanifold_components(M, p)) + e = [(t[1], check_size(t[2:end]...)) for t in ts] + errors = filter((x) -> !(x[2] === nothing), e) + cerr = [ComponentManifoldError(er...) for er in errors] + (length(errors) > 1) && return CompositeManifoldError(cerr) + (length(errors) == 1) && return cerr[1] + return nothing +end + +function check_size(M::ProductManifold, p, X) + try + submanifold_components(M, X) + catch e + return DomainError("Vector $X does not support submanifold_components") + end + ts = ziptuples( + Tuple(1:length(M.manifolds)), + M.manifolds, + submanifold_components(M, p), + submanifold_components(M, X), + ) + e = [(t[1], check_size(t[2:end]...)) for t in ts] + errors = filter(x -> !(x[2] === nothing), e) + cerr = [ComponentManifoldError(er...) for er in errors] + (length(errors) > 1) && return CompositeManifoldError(cerr) + (length(errors) == 1) && return cerr[1] + return nothing +end +""" + check_vector(M::ProductManifold, p, X; kwargs... ) + +Check whether `X` is a tangent vector to `p` on the [`ProductManifold`](@ref) +`M`, i.e. all projections to base manifolds must be respective tangent vectors. +If `X` is not a tangent vector to `p` on `M` a [`CompositeManifoldError`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/functions.html#ManifoldsBase.CompositeManifoldError).consisting +of all error messages of the components, for which the tests fail is returned. + +The tolerance for the last test can be set using the `kwargs...`. +""" +function check_vector(M::ProductManifold, p, X; kwargs...) + try + submanifold_components(M, X) + catch e + return DomainError("Vector $X does not support submanifold_components") + end + ts = ziptuples( + Tuple(1:length(M.manifolds)), + M.manifolds, + submanifold_components(M, p), + submanifold_components(M, X), + ) + e = [(t[1], check_vector(t[2:end]...; kwargs...)) for t in ts] + errors = filter(x -> !(x[2] === nothing), e) + cerr = [ComponentManifoldError(er...) for er in errors] + (length(errors) > 1) && return CompositeManifoldError(cerr) + (length(errors) == 1) && return cerr[1] + return nothing +end + +@doc raw""" + ×(M, N) + cross(M, N) + cross(M1, M2, M3,...) + +Return the [`ProductManifold`](@ref) For two `AbstractManifold`s `M` and `N`, +where for the case that one of them is a [`ProductManifold`](@ref) itself, +the other is either prepended (if `N` is a product) or appenden (if `M`) is. +If both are product manifold, they are combined into one product manifold, +keeping the order. + +For the case that more than one is a product manifold of these is build with the +same approach as above +""" +cross(::AbstractManifold...) +LinearAlgebra.cross(M::AbstractManifold, N::AbstractManifold) = ProductManifold(M, N) +function LinearAlgebra.cross(M::ProductManifold, N::AbstractManifold) + return ProductManifold(M.manifolds..., N) +end +function LinearAlgebra.cross(M::AbstractManifold, N::ProductManifold) + return ProductManifold(M, N.manifolds...) +end +function LinearAlgebra.cross(M::ProductManifold, N::ProductManifold) + return ProductManifold(M.manifolds..., N.manifolds...) +end + + +@doc raw""" + ×(m, n) + cross(m, n) + cross(m1, m2, m3,...) + +Return the [`ProductRetraction`](@ref) For two or more [`AbstractRetractionMethod`](@ref)s, +where for the case that one of them is a [`ProductRetraction`](@ref) itself, +the other is either prepended (if `m` is a product) or appenden (if `n`) is. +If both [`ProductRetraction`](@ref)s, they are combined into one keeping the order. +""" +cross(::AbstractRetractionMethod...) +function LinearAlgebra.cross(m::AbstractRetractionMethod, n::AbstractRetractionMethod) + return ProductRetraction(m, n) +end +function LinearAlgebra.cross(m::ProductRetraction, n::AbstractRetractionMethod) + return ProductRetraction(m.retractions..., n) +end +function LinearAlgebra.cross(m::AbstractRetractionMethod, n::ProductRetraction) + return ProductRetraction(m, n.retractions...) +end +function LinearAlgebra.cross(m::ProductRetraction, n::ProductRetraction) + return ProductRetraction(m.retractions..., n.retractions...) +end + +@doc raw""" + ×(m, n) + cross(m, n) + cross(m1, m2, m3,...) + +Return the [`InverseProductRetraction`](@ref) For two or more [`AbstractInverseRetractionMethod`](@ref)s, +where for the case that one of them is a [`InverseProductRetraction`](@ref) itself, +the other is either prepended (if `r` is a product) or appenden (if `s`) is. +If both [`InverseProductRetraction`](@ref)s, they are combined into one keeping the order. +""" +cross(::AbstractInverseRetractionMethod...) +function LinearAlgebra.cross( + m::AbstractInverseRetractionMethod, + n::AbstractInverseRetractionMethod, +) + return InverseProductRetraction(m, n) +end +function LinearAlgebra.cross( + m::InverseProductRetraction, + n::AbstractInverseRetractionMethod, +) + return InverseProductRetraction(m.inverse_retractions..., n) +end +function LinearAlgebra.cross( + m::AbstractInverseRetractionMethod, + n::InverseProductRetraction, +) + return InverseProductRetraction(m, n.inverse_retractions...) +end +function LinearAlgebra.cross(m::InverseProductRetraction, n::InverseProductRetraction) + return InverseProductRetraction(m.inverse_retractions..., n.inverse_retractions...) +end + + +@doc raw""" + ×(m, n) + cross(m, n) + cross(m1, m2, m3,...) + +Return the [`ProductVectorTransport`](@ref) For two or more [`AbstractVectorTransportMethod`](@ref)s, +where for the case that one of them is a [`ProductVectorTransport`](@ref) itself, +the other is either prepended (if `r` is a product) or appenden (if `s`) is. +If both [`ProductVectorTransport`](@ref)s, they are combined into one keeping the order. +""" +cross(::AbstractVectorTransportMethod...) +function LinearAlgebra.cross( + m::AbstractVectorTransportMethod, + n::AbstractVectorTransportMethod, +) + return ProductVectorTransport(m, n) +end +function LinearAlgebra.cross(m::ProductVectorTransport, n::AbstractVectorTransportMethod) + return ProductVectorTransport(m.methods..., n) +end +function LinearAlgebra.cross(m::AbstractVectorTransportMethod, n::ProductVectorTransport) + return ProductVectorTransport(m, n.methods...) +end +function LinearAlgebra.cross(m::ProductVectorTransport, n::ProductVectorTransport) + return ProductVectorTransport(m.methods..., n.methods...) +end + +function default_retraction_method(M::ProductManifold) + return ProductRetraction(map(default_retraction_method, M.manifolds)...) +end + +function default_inverse_retraction_method(M::ProductManifold) + return InverseProductRetraction(map(default_inverse_retraction_method, M.manifolds)...) +end + +function default_vector_transport_method(M::ProductManifold) + return ProductVectorTransport(map(default_vector_transport_method, M.manifolds)...) +end + +@doc raw""" + distance(M::ProductManifold, p, q) + +Compute the distance between two points `p` and `q` on the [`ProductManifold`](@ref) `M`, which is +the 2-norm of the elementwise distances on the internal manifolds that build `M`. +""" +function distance(M::ProductManifold, p, q) + return sqrt( + sum( + map( + distance, + M.manifolds, + submanifold_components(M, p), + submanifold_components(M, q), + ) .^ 2, + ), + ) +end + +@doc raw""" + exp(M::ProductManifold, p, X) + +compute the exponential map from `p` in the direction of `X` on the [`ProductManifold`](@ref) `M`, +which is the elementwise exponential map on the internal manifolds that build `M`. +""" +exp(::ProductManifold, ::Any...) + +function exp!(M::ProductManifold, q, p, X) + map( + exp!, + M.manifolds, + submanifold_components(M, q), + submanifold_components(M, p), + submanifold_components(M, X), + ) + return q +end +function exp!(M::ProductManifold, q, p, X, t::Number) + map( + (N, qc, pc, Xc) -> exp!(N, qc, pc, Xc, t), + M.manifolds, + submanifold_components(M, q), + submanifold_components(M, p), + submanifold_components(M, X), + ) + return q +end + +function get_basis(M::ProductManifold, p, B::AbstractBasis) + parts = map(t -> get_basis(t..., B), ziptuples(M.manifolds, submanifold_components(p))) + return CachedBasis(B, ProductBasisData(parts)) +end +function get_basis(M::ProductManifold, p, B::CachedBasis) + return invoke(get_basis, Tuple{AbstractManifold,Any,CachedBasis}, M, p, B) +end +function get_basis(M::ProductManifold, p, B::DiagonalizingOrthonormalBasis) + vs = map( + ziptuples( + M.manifolds, + submanifold_components(p), + submanifold_components(B.frame_direction), + ), + ) do t + return get_basis(t[1], t[2], DiagonalizingOrthonormalBasis(t[3])) + end + return CachedBasis(B, ProductBasisData(vs)) +end + +""" + get_component(M::ProductManifold, p, i) + +Get the `i`th component of a point `p` on a [`ProductManifold`](@ref) `M`. +""" +function get_component(M::ProductManifold, p, i) + return submanifold_component(M, p, i) +end + +function get_coordinates(M::ProductManifold, p, X, B::AbstractBasis) + reps = map( + t -> get_coordinates(t..., B), + ziptuples(M.manifolds, submanifold_components(M, p), submanifold_components(M, X)), + ) + return vcat(reps...) +end +function get_coordinates( + M::ProductManifold, + p, + X, + B::CachedBasis{𝔽,<:AbstractBasis{𝔽},<:ProductBasisData}, +) where {𝔽} + reps = map( + get_coordinates, + M.manifolds, + submanifold_components(M, p), + submanifold_components(M, X), + B.data.parts, + ) + return vcat(reps...) +end + +function get_coordinates!(M::ProductManifold, Xⁱ, p, X, B::AbstractBasis) + dim = manifold_dimension(M) + @assert length(Xⁱ) == dim + i = one(dim) + ts = ziptuples(M.manifolds, submanifold_components(M, p), submanifold_components(M, X)) + for t in ts + SM = first(t) + dim = manifold_dimension(SM) + tXⁱ = @inbounds view(Xⁱ, i:(i + dim - 1)) + get_coordinates!(SM, tXⁱ, Base.tail(t)..., B) + i += dim + end + return Xⁱ +end +function get_coordinates!( + M::ProductManifold, + Xⁱ, + p, + X, + B::CachedBasis{𝔽,<:AbstractBasis{𝔽},<:ProductBasisData}, +) where {𝔽} + dim = manifold_dimension(M) + @assert length(Xⁱ) == dim + i = one(dim) + ts = ziptuples( + M.manifolds, + submanifold_components(M, p), + submanifold_components(M, X), + B.data.parts, + ) + for t in ts + SM = first(t) + dim = manifold_dimension(SM) + tXⁱ = @inbounds view(Xⁱ, i:(i + dim - 1)) + get_coordinates!(SM, tXⁱ, Base.tail(t)...) + i += dim + end + return Xⁱ +end + +function _get_dim_ranges(dims::NTuple{N,Any}) where {N} + dims_acc = accumulate(+, (1, dims...)) + return ntuple(i -> (dims_acc[i]:(dims_acc[i] + dims[i] - 1)), Val(N)) +end + +function get_vector!(M::ProductManifold, X, p, Xⁱ, B::AbstractBasis) + dims = map(manifold_dimension, M.manifolds) + @assert length(Xⁱ) == sum(dims) + dim_ranges = _get_dim_ranges(dims) + tXⁱ = map(dr -> (@inbounds view(Xⁱ, dr)), dim_ranges) + ts = ziptuples( + M.manifolds, + submanifold_components(M, X), + submanifold_components(M, p), + tXⁱ, + ) + map(ts) do t + return get_vector!(t..., B) + end + return X +end +function get_vector!( + M::ProductManifold, + X, + p, + Xⁱ, + B::CachedBasis{𝔽,<:AbstractBasis{𝔽},<:ProductBasisData}, +) where {𝔽} + dims = map(manifold_dimension, M.manifolds) + @assert length(Xⁱ) == sum(dims) + dim_ranges = _get_dim_ranges(dims) + tXⁱ = map(dr -> (@inbounds view(Xⁱ, dr)), dim_ranges) + ts = ziptuples( + M.manifolds, + submanifold_components(M, X), + submanifold_components(M, p), + tXⁱ, + B.data.parts, + ) + map(ts) do t + return get_vector!(t...) + end + return X +end + +@doc raw""" + injectivity_radius(M::ProductManifold) + injectivity_radius(M::ProductManifold, x) + +Compute the injectivity radius on the [`ProductManifold`](@ref), which is the +minimum of the factor manifolds. +""" +injectivity_radius(::ProductManifold, ::Any...) +function injectivity_radius(M::ProductManifold, p) + return min(map(injectivity_radius, M.manifolds, submanifold_components(M, p))...) +end +function injectivity_radius(M::ProductManifold, p, m::AbstractRetractionMethod) + return min( + map( + (lM, lp) -> injectivity_radius(lM, lp, m), + M.manifolds, + submanifold_components(M, p), + )..., + ) +end +function injectivity_radius(M::ProductManifold, p, m::ProductRetraction) + return min( + map( + (lM, lp, lm) -> injectivity_radius(lM, lp, lm), + M.manifolds, + submanifold_components(M, p), + m.retractions, + )..., + ) +end +injectivity_radius(M::ProductManifold) = min(map(injectivity_radius, M.manifolds)...) +function injectivity_radius(M::ProductManifold, m::AbstractRetractionMethod) + return min(map(manif -> injectivity_radius(manif, m), M.manifolds)...) +end +function injectivity_radius(M::ProductManifold, m::ProductRetraction) + return min(map((lM, lm) -> injectivity_radius(lM, lm), M.manifolds, m.retractions)...) +end + +@doc raw""" + inner(M::ProductManifold, p, X, Y) + +compute the inner product of two tangent vectors `X`, `Y` from the tangent space +at `p` on the [`ProductManifold`](@ref) `M`, which is just the sum of the +internal manifolds that build `M`. +""" +function inner(M::ProductManifold, p, X, Y) + subproducts = map( + inner, + M.manifolds, + submanifold_components(M, p), + submanifold_components(M, X), + submanifold_components(M, Y), + ) + return sum(subproducts) +end + +@doc raw""" + inverse_retract(M::ProductManifold, p, q, m::InverseProductRetraction) + +Compute the inverse retraction from `p` with respect to `q` on the [`ProductManifold`](@ref) +`M` using an [`InverseProductRetraction`](@ref), which by default encapsulates a inverse +retraction for each manifold of the product. Then this method is performed elementwise, +so the encapsulated inverse retraction methods have to be available per factor. +""" +inverse_retract(::ProductManifold, ::Any, ::Any, ::Any, ::InverseProductRetraction) + +@doc raw""" + inverse_retract(M::ProductManifold, p, q, m::AbstractInverseRetractionMethod) + +Compute the inverse retraction from `p` with respect to `q` on the [`ProductManifold`](@ref) +`M` using an [`AbstractInverseRetractionMethod`](@ref), which is used on each manifold of +the product. +""" +inverse_retract(::ProductManifold, ::Any, ::Any, ::Any, ::AbstractInverseRetractionMethod) + +function inverse_retract!(M::ProductManifold, Y, p, q, method::InverseProductRetraction) + map( + (iM, iY, ip, iq, im) -> inverse_retract!(iM, iY, ip, iq, im), + M.manifolds, + submanifold_components(M, Y), + submanifold_components(M, p), + submanifold_components(M, q), + method.inverse_retractions, + ) + return Y +end +function inverse_retract!( + M::ProductManifold, + Y, + p, + q, + method::IRM, +) where {IRM<:AbstractInverseRetractionMethod} + map( + (iM, iY, ip, iq) -> inverse_retract!(iM, iY, ip, iq, method), + M.manifolds, + submanifold_components(M, Y), + submanifold_components(M, p), + submanifold_components(M, q), + ) + return Y +end +function _isapprox(M::ProductManifold, p, q; kwargs...) + return all( + t -> isapprox(t...; kwargs...), + ziptuples(M.manifolds, submanifold_components(M, p), submanifold_components(M, q)), + ) +end +function _isapprox(M::ProductManifold, p, X, Y; kwargs...) + return all( + t -> isapprox(t...; kwargs...), + ziptuples( + M.manifolds, + submanifold_components(M, p), + submanifold_components(M, X), + submanifold_components(M, Y), + ), + ) +end + +""" + is_flat(::ProductManifold) + +Return true if and only if all component manifolds of [`ProductManifold`](@ref) `M` are flat. +""" +function is_flat(M::ProductManifold) + return all(is_flat, M.manifolds) +end + +@doc raw""" + log(M::ProductManifold, p, q) + +Compute the logarithmic map from `p` to `q` on the [`ProductManifold`](@ref) `M`, +which can be computed using the logarithmic maps of the manifolds elementwise. +""" +log(::ProductManifold, ::Any...) + +function log!(M::ProductManifold, X, p, q) + map( + log!, + M.manifolds, + submanifold_components(M, X), + submanifold_components(M, p), + submanifold_components(M, q), + ) + return X +end + +@doc raw""" + manifold_dimension(M::ProductManifold) + +Return the manifold dimension of the [`ProductManifold`](@ref), which is the sum of the +manifold dimensions the product is made of. +""" +manifold_dimension(M::ProductManifold) = mapreduce(manifold_dimension, +, M.manifolds) + +function mid_point!(M::ProductManifold, q, p1, p2) + map( + mid_point!, + M.manifolds, + submanifold_components(M, q), + submanifold_components(M, p1), + submanifold_components(M, p2), + ) + return q +end + +@doc raw""" + norm(M::ProductManifold, p, X) + +Compute the norm of `X` from the tangent space of `p` on the [`ProductManifold`](@ref), +i.e. from the element wise norms the 2-norm is computed. +""" +function LinearAlgebra.norm(M::ProductManifold, p, X) + norms_squared = ( + map( + norm, + M.manifolds, + submanifold_components(M, p), + submanifold_components(M, X), + ) .^ 2 + ) + return sqrt(sum(norms_squared)) +end + +""" + number_of_components(M::ProductManifold{<:NTuple{N,Any}}) where {N} + +Calculate the number of manifolds multiplied in the given [`ProductManifold`](@ref) `M`. +""" +number_of_components(::ProductManifold{𝔽,<:NTuple{N,Any}}) where {𝔽,N} = N + +function parallel_transport_direction!(M::ProductManifold, Y, p, X, d) + map( + parallel_transport_direction!, + M.manifolds, + submanifold_components(M, Y), + submanifold_components(M, p), + submanifold_components(M, X), + submanifold_components(M, d), + ) + return Y +end + +function parallel_transport_to!(M::ProductManifold, Y, p, X, q) + map( + parallel_transport_to!, + M.manifolds, + submanifold_components(M, Y), + submanifold_components(M, p), + submanifold_components(M, X), + submanifold_components(M, q), + ) + return Y +end + +function project!(M::ProductManifold, q, p) + map(project!, M.manifolds, submanifold_components(M, q), submanifold_components(M, p)) + return q +end + +function project!(M::ProductManifold, Y, p, X) + map( + project!, + M.manifolds, + submanifold_components(M, Y), + submanifold_components(M, p), + submanifold_components(M, X), + ) + return Y +end + +function Random.rand!( + M::ProductManifold, + pX; + vector_at = nothing, + parts_kwargs = map(_ -> (;), M.manifolds), +) + return rand!( + Random.default_rng(), + M, + pX; + vector_at = vector_at, + parts_kwargs = parts_kwargs, + ) +end +function Random.rand!( + rng::AbstractRNG, + M::ProductManifold, + pX; + vector_at = nothing, + parts_kwargs = map(_ -> (;), M.manifolds), +) + if vector_at === nothing + map( + (N, q, kwargs) -> rand!(rng, N, q; kwargs...), + M.manifolds, + submanifold_components(M, pX), + parts_kwargs, + ) + else + map( + (N, X, p, kwargs) -> rand!(rng, N, X; vector_at = p, kwargs...), + M.manifolds, + submanifold_components(M, pX), + submanifold_components(M, vector_at), + parts_kwargs, + ) + end + return pX +end + +@doc raw""" + retract(M::ProductManifold, p, X, m::ProductRetraction) + +Compute the retraction from `p` with tangent vector `X` on the [`ProductManifold`](@ref) `M` +using an [`ProductRetraction`](@ref), which by default encapsulates retractions of the +base manifolds. Then this method is performed elementwise, so the encapsulated retractions +method has to be one that is available on the manifolds. +""" +retract(::ProductManifold, ::Any, ::Any, ::ProductRetraction) + +@doc raw""" + retract(M::ProductManifold, p, X, m::AbstractRetractionMethod) + +Compute the retraction from `p` with tangent vector `X` on the [`ProductManifold`](@ref) `M` +using the [`AbstractRetractionMethod`](@ref) `m` on every manifold. +""" +retract(::ProductManifold, ::Any, ::Any, ::AbstractRetractionMethod) + +function retract!(M::ProductManifold, q, p, X, t::Number, method::ProductRetraction) + map( + (N, qc, pc, Xc, rm) -> retract!(N, qc, pc, Xc, t, rm), + M.manifolds, + submanifold_components(M, q), + submanifold_components(M, p), + submanifold_components(M, X), + method.retractions, + ) + return q +end +function retract!( + M::ProductManifold, + q, + p, + X, + t::Number, + method::RTM, +) where {RTM<:AbstractRetractionMethod} + map( + (N, qc, pc, Xc) -> retract!(N, qc, pc, Xc, t, method), + M.manifolds, + submanifold_components(M, q), + submanifold_components(M, p), + submanifold_components(M, X), + ) + return q +end + +function representation_size(M::ProductManifold) + return (mapreduce(m -> prod(representation_size(m)), +, M.manifolds),) +end + +@doc raw""" + riemann_tensor(M::ProductManifold, p, X, Y, Z) + +Compute the Riemann tensor at point from `p` with tangent vectors `X`, `Y` and `Z` on +the [`ProductManifold`](@ref) `M`. +""" +riemann_tensor(M::ProductManifold, p, X, Y, X) + +function riemann_tensor!(M::ProductManifold, Xresult, p, X, Y, Z) + map( + riemann_tensor!, + M.manifolds, + submanifold_components(M, Xresult), + submanifold_components(M, p), + submanifold_components(M, X), + submanifold_components(M, Y), + submanifold_components(M, Z), + ) + return Xresult +end + +""" + select_from_tuple(t::NTuple{N, Any}, positions::Val{P}) + +Selects elements of tuple `t` at positions specified by the second argument. +For example `select_from_tuple(("a", "b", "c"), Val((3, 1, 1)))` returns +`("c", "a", "a")`. +""" +@generated function select_from_tuple(t::NTuple{N,Any}, positions::Val{P}) where {N,P} + for k in P + (k < 0 || k > N) && error("positions must be between 1 and $N") + end + return Expr(:tuple, [Expr(:ref, :t, k) for k in P]...) +end + +""" + set_component!(M::ProductManifold, q, p, i) + +Set the `i`th component of a point `q` on a [`ProductManifold`](@ref) `M` to `p`, where `p` is a point on the [`AbstractManifold`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/types.html#ManifoldsBase.AbstractManifold) this factor of the product manifold consists of. +""" +function set_component!(M::ProductManifold, q, p, i) + return copyto!(submanifold_component(M, q, i), p) +end + +function _show_submanifold(io::IO, M::AbstractManifold; pre = "") + sx = sprint(show, "text/plain", M, context = io, sizehint = 0) + if occursin('\n', sx) + sx = sprint(show, M, context = io, sizehint = 0) + end + sx = replace(sx, '\n' => "\n$(pre)") + print(io, pre, sx) + return nothing +end + +function _show_submanifold_range(io::IO, Ms, range; pre = "") + for i in range + M = Ms[i] + print(io, '\n') + _show_submanifold(io, M; pre = pre) + end + return nothing +end + +function _show_product_manifold_no_header(io::IO, M) + n = length(M.manifolds) + sz = displaysize(io) + screen_height, screen_width = sz[1] - 4, sz[2] + half_height = div(screen_height, 2) + inds = 1:n + pre = " " + if n > screen_height + inds = [1:half_height; (n - div(screen_height - 1, 2) + 1):n] + end + if n ≤ screen_height + _show_submanifold_range(io, M.manifolds, 1:n; pre = pre) + else + _show_submanifold_range(io, M.manifolds, 1:half_height; pre = pre) + print(io, "\n$(pre)⋮") + _show_submanifold_range( + io, + M.manifolds, + (n - div(screen_height - 1, 2) + 1):n; + pre = pre, + ) + end + return nothing +end + +function Base.show(io::IO, ::MIME"text/plain", M::ProductManifold) + n = length(M.manifolds) + print(io, "ProductManifold with $(n) submanifold$(n == 1 ? "" : "s"):") + return _show_product_manifold_no_header(io, M) +end + +function Base.show(io::IO, M::ProductManifold) + return print(io, "ProductManifold(", join(M.manifolds, ", "), ")") +end + +function Base.show(io::IO, m::ProductRetraction) + return print(io, "ProductRetraction(", join(m.retractions, ", "), ")") +end + +function Base.show(io::IO, m::InverseProductRetraction) + return print(io, "InverseProductRetraction(", join(m.inverse_retractions, ", "), ")") +end + +function Base.show(io::IO, m::ProductVectorTransport) + return print(io, "ProductVectorTransport(", join(m.methods, ", "), ")") +end + + +function Base.show( + io::IO, + mime::MIME"text/plain", + B::CachedBasis{𝔽,T,D}, +) where {𝔽,T<:AbstractBasis{𝔽},D<:ProductBasisData} + println(io, "$(T) for a product manifold") + for (i, cb) in enumerate(B.data.parts) + println(io, "Basis for component $i:") + show(io, mime, cb) + println(io) + end + return nothing +end + +""" + submanifold(M::ProductManifold, i::Integer) + +Extract the `i`th factor of the product manifold `M`. +""" +submanifold(M::ProductManifold, i::Integer) = M.manifolds[i] + +""" + submanifold(M::ProductManifold, i::Val) + submanifold(M::ProductManifold, i::AbstractVector) + +Extract the factor of the product manifold `M` indicated by indices in `i`. +For example, for `i` equal to `Val((1, 3))` the product manifold constructed +from the first and the third factor is returned. + +The version with `AbstractVector` is not type-stable, for better preformance use `Val`. +""" +function submanifold(M::ProductManifold, i::Val) + return ProductManifold(select_from_tuple(M.manifolds, i)...) +end +submanifold(M::ProductManifold, i::AbstractVector) = submanifold(M, Val(tuple(i...))) + +function vector_transport_direction!( + M::ProductManifold, + Y, + p, + X, + d, + m::ProductVectorTransport, +) + map( + vector_transport_direction!, + M.manifolds, + submanifold_components(M, Y), + submanifold_components(M, p), + submanifold_components(M, X), + submanifold_components(M, d), + m.methods, + ) + return Y +end +function vector_transport_direction!( + M::ProductManifold, + Y, + p, + X, + d, + m::VTM, +) where {VTM<:AbstractVectorTransportMethod} + map( + (iM, iY, ip, iX, id) -> vector_transport_direction!(iM, iY, ip, iX, id, m), + M.manifolds, + submanifold_components(M, Y), + submanifold_components(M, p), + submanifold_components(M, X), + submanifold_components(M, d), + ) + return Y +end + +@doc raw""" + vector_transport_to(M::ProductManifold, p, X, q, m::ProductVectorTransport) + +Compute the vector transport the tangent vector `X` at `p` to `q` on the +[`ProductManifold`](@ref) `M` using a [`ProductVectorTransport`](@ref) `m`. +""" +vector_transport_to(::ProductManifold, ::Any, ::Any, ::Any, ::ProductVectorTransport) + +@doc raw""" + vector_transport_to(M::ProductManifold, p, X, q, m::AbstractVectorTransportMethod) + +Compute the vector transport the tangent vector `X` at `p` to `q` on the +[`ProductManifold`](@ref) `M` using an [`AbstractVectorTransportMethod`](@ref) `m` +on each manifold. +""" +vector_transport_to(::ProductManifold, ::Any, ::Any, ::Any, ::AbstractVectorTransportMethod) + + +function vector_transport_to!(M::ProductManifold, Y, p, X, q, m::ProductVectorTransport) + map( + vector_transport_to!, + M.manifolds, + submanifold_components(M, Y), + submanifold_components(M, p), + submanifold_components(M, X), + submanifold_components(M, q), + m.methods, + ) + return Y +end +function vector_transport_to!( + M::ProductManifold, + Y, + p, + X, + q, + m::AbstractVectorTransportMethod, +) + map( + (iM, iY, ip, iX, iq) -> vector_transport_to!(iM, iY, ip, iX, iq, m), + M.manifolds, + submanifold_components(M, Y), + submanifold_components(M, p), + submanifold_components(M, X), + submanifold_components(M, q), + ), + return Y +end + +@doc raw""" + Y = Weingarten(M::ProductManifold, p, X, V) + Weingarten!(M::ProductManifold, Y, p, X, V) + +Since the metric decouples, also the computation of the Weingarten map +``\mathcal W_p`` can be computed elementwise on the single elements of the [`ProductManifold`](@ref) `M`. +""" +Weingarten(::ProductManifold, p, X, V) + +function Weingarten!(M::ProductManifold, Y, p, X, V) + map( + Weingarten!, + M.manifolds, + submanifold_components(M, Y), + submanifold_components(M, p), + submanifold_components(M, X), + submanifold_components(M, V), + ) + return Y +end + +function zero_vector!(M::ProductManifold, X, p) + map( + zero_vector!, + M.manifolds, + submanifold_components(M, X), + submanifold_components(M, p), + ) + return X +end + + + + +@doc raw""" + submanifold_component(M::AbstractManifold, p, i::Integer) + submanifold_component(M::AbstractManifold, p, ::Val{i}) where {i} + submanifold_component(p, i::Integer) + submanifold_component(p, ::Val{i}) where {i} + +Project the product array `p` on `M` to its `i`th component. A new array is returned. +""" +submanifold_component(::Any...) +@inline function submanifold_component(M::AbstractManifold, p, i::Integer) + return submanifold_component(M, p, Val(i)) +end +@inline submanifold_component(M::AbstractManifold, p, i::Val) = submanifold_component(p, i) +@inline submanifold_component(p, i::Integer) = submanifold_component(p, Val(i)) + +@doc raw""" + submanifold_components(M::AbstractManifold, p) + submanifold_components(p) + +Get the projected components of `p` on the submanifolds of `M`. The components are returned in a Tuple. +""" +submanifold_components(::Any...) +@inline submanifold_components(::AbstractManifold, p) = submanifold_components(p) + +""" + ziptuples(a, b[, c[, d[, e]]]) + +Zips tuples `a`, `b`, and remaining in a fast, type-stable way. If they have different +lengths, the result is trimmed to the length of the shorter tuple. +""" +@generated function ziptuples(a::NTuple{N,Any}, b::NTuple{M,Any}) where {N,M} + ex = Expr(:tuple) + for i in 1:min(N, M) + push!(ex.args, :((a[$i], b[$i]))) + end + return ex +end +@generated function ziptuples( + a::NTuple{N,Any}, + b::NTuple{M,Any}, + c::NTuple{L,Any}, +) where {N,M,L} + ex = Expr(:tuple) + for i in 1:min(N, M, L) + push!(ex.args, :((a[$i], b[$i], c[$i]))) + end + return ex +end +@generated function ziptuples( + a::NTuple{N,Any}, + b::NTuple{M,Any}, + c::NTuple{L,Any}, + d::NTuple{K,Any}, +) where {N,M,L,K} + ex = Expr(:tuple) + for i in 1:min(N, M, L, K) + push!(ex.args, :((a[$i], b[$i], c[$i], d[$i]))) + end + return ex +end +@generated function ziptuples( + a::NTuple{N,Any}, + b::NTuple{M,Any}, + c::NTuple{L,Any}, + d::NTuple{K,Any}, + e::NTuple{J,Any}, +) where {N,M,L,K,J} + ex = Expr(:tuple) + for i in 1:min(N, M, L, K, J) + push!(ex.args, :((a[$i], b[$i], c[$i], d[$i], e[$i]))) + end + return ex +end diff --git a/src/TangentSpace.jl b/src/TangentSpace.jl new file mode 100644 index 00000000..6e568b25 --- /dev/null +++ b/src/TangentSpace.jl @@ -0,0 +1,254 @@ + +@doc raw""" + TangentSpace{𝔽,M} = Fiber{𝔽,TangentSpaceType,M} where {𝔽,M<:AbstractManifold{𝔽}} + +A manifold for the tangent space ``T_p\mathcal M`` at a point ``p\in\mathcal M``. +This is modelled as an alias for [`VectorSpaceFiber`](@ref) corresponding to +[`TangentSpaceType`](@ref). + +# Constructor + + TangentSpace(M::AbstractManifold, p) + +Return the manifold (vector space) representing the tangent space ``T_p\mathcal M`` +at point `p`, ``p\in\mathcal M``. +""" +const TangentSpace{𝔽,M} = Fiber{𝔽,TangentSpaceType,M} where {𝔽,M<:AbstractManifold{𝔽}} + +TangentSpace(M::AbstractManifold, p) = Fiber(M, p, TangentSpaceType()) + +@doc raw""" + CotangentSpace{𝔽,M} = Fiber{𝔽,CotangentSpaceType,M} where {𝔽,M<:AbstractManifold{𝔽}} + +A manifold for the Cotangent space ``T^*_p\mathcal M`` at a point ``p\in\mathcal M``. +This is modelled as an alias for [`VectorSpaceFiber`](@ref) corresponding to +[`CotangentSpaceType`](@ref). + +# Constructor + + CotangentSpace(M::AbstractManifold, p) + +Return the manifold (vector space) representing the cotangent space ``T^*_p\mathcal M`` +at point `p`, ``p\in\mathcal M``. +""" +const CotangentSpace{𝔽,M} = Fiber{𝔽,CotangentSpaceType,M} where {𝔽,M<:AbstractManifold{𝔽}} + +CotangentSpace(M::AbstractManifold, p) = Fiber(M, p, CotangentSpaceType()) + + +function allocate_result(M::TangentSpace, ::typeof(rand)) + return zero_vector(M.manifold, M.point) +end + +""" + distance(M::TangentSpace, X, Y) + +Distance between vectors `X` and `Y` from the [`TangentSpace`](@ref) `TpM`. +It is calculated as the [`norm`](@ref) (induced by the metric on `TpM`) of their difference. +""" +function distance(TpM::TangentSpace, X, Y) + return norm(TpM.manifold, TpM.point, Y - X) +end + +function embed!(TpM::TangentSpace, Y, X) + return embed!(TpM.manifold, Y, TpM.point, X) +end +function embed!(TpM::TangentSpace, W, X, V) + return embed!(TpM.manifold, W, TpM.point, V) +end + +@doc raw""" + exp(TpM::TangentSpace, X, V) + +Exponential map of tangent vectors `X` from `TpM` and a direction `V`, +which is also from the [`TangentSpace`](@ref) `TpM` since we identify the tangent space of `TpM` with `TpM`. +The exponential map then simplifies to the sum `X+V`. +""" +exp(::TangentSpace, ::Any, ::Any) + +function exp!(TpM::TangentSpace, Y, X, V) + copyto!(TpM.manifold, Y, TpM.point, X + V) + return Y +end + +fiber_dimension(M::AbstractManifold, ::CotangentSpaceType) = manifold_dimension(M) +fiber_dimension(M::AbstractManifold, ::TangentSpaceType) = manifold_dimension(M) + +function get_basis(TpM::TangentSpace, X, B::CachedBasis) + return invoke( + get_basis, + Tuple{AbstractManifold,Any,CachedBasis}, + TpM.manifold, + TpM.point, + B, + ) +end +function get_basis(TpM::TangentSpace, X, B::AbstractBasis{<:Any,TangentSpaceType}) + return get_basis(TpM.manifold, TpM.point, B) +end + +function get_coordinates(TpM::TangentSpace, X, V, B::AbstractBasis) + return get_coordinates(TpM.manifold, TpM.point, V, B) +end + +function get_coordinates!(TpM::TangentSpace, c, X, V, B::AbstractBasis) + return get_coordinates!(TpM.manifold, c, TpM.point, V, B) +end + +function get_vector(TpM::TangentSpace, X, c, B::AbstractBasis) + return get_vector(TpM.manifold, TpM.point, c, B) +end + +function get_vector!(TpM::TangentSpace, V, X, c, B::AbstractBasis) + return get_vector!(TpM.manifold, V, TpM.point, c, B) +end + +function get_vectors(TpM::TangentSpace, X, B::CachedBasis) + return get_vectors(TpM.manifold, TpM.point, B) +end + +@doc raw""" + injectivity_radius(TpM::TangentSpace) + +Return the injectivity radius on the [`TangentSpace`](@ref) `TpM`, which is $∞$. +""" +injectivity_radius(::TangentSpace) = Inf + +@doc raw""" + inner(M::TangentSpace, X, V, W) + +For any ``X ∈ T_p\mathcal M`` we identify the tangent space ``T_X(T_p\mathcal M)`` +with ``T_p\mathcal M`` again. Hence an inner product of ``V,W`` is just the inner product of +the tangent space itself. ``⟨V,W⟩_X = ⟨V,W⟩_p``. +""" +function inner(TpM::TangentSpace, X, V, W) + return inner(TpM.manifold, TpM.point, V, W) +end + +""" + is_flat(::TangentSpace) + +The [`TangentSpace`](@ref) is a flat manifold, so this returns `true`. +""" +is_flat(::TangentSpace) = true + +function _isapprox(TpM::TangentSpace, X, Y; kwargs...) + return isapprox(TpM.manifold, TpM.point, X, Y; kwargs...) +end + +function _isapprox(TpM::TangentSpace, X, V, W; kwargs...) + return isapprox(TpM.manifold, TpM.point, V, W; kwargs...) +end + +""" + log(TpM::TangentSpace, X, Y) + +Logarithmic map on the [`TangentSpace`](@ref) `TpM`, calculated as the difference of tangent +vectors `q` and `p` from `TpM`. +""" +log(::TangentSpace, ::Any...) +function log!(TpM::TangentSpace, V, X, Y) + copyto!(TpM, V, TpM.point, Y - X) + return V +end + +@doc raw""" + manifold_dimension(TpM::TangentSpace) + +Return the dimension of the [`TangentSpace`](@ref) ``T_p\mathcal M`` at ``p∈\mathcal M``, +which is the same as the dimension of the manifold ``\mathcal M``. +""" +function manifold_dimension(TpM::TangentSpace) + return manifold_dimension(TpM.manifold) +end + +@doc raw""" + parallel_transport_to(::TangentSpace, X, V, Y) + +Transport the tangent vector ``Z ∈ T_X(T_p\mathcal M)`` from `X` to `Y`. +Since we identify ``T_X(T_p\mathcal M) = T_p\mathcal M`` and the tangent space is a vector space, +parallel transport simplifies to the identity, so this function yields ``V`` as a result. +""" +parallel_transport_to(TpM::TangentSpace, X, V, Y) + +function parallel_transport_to!(TpM::TangentSpace, W, X, V, Y) + return copyto!(TpM.manifold, W, TpM.point, V) +end + +@doc raw""" + project(TpM::TangentSpace, X) + +Project the point `X` from embedding of the [`TangentSpace`](@ref) `TpM` onto `TpM`. +""" +project(::TangentSpace, ::Any) + +function project!(TpM::TangentSpace, Y, X) + return project!(TpM.manifold, Y, TpM.point, X) +end + +@doc raw""" + project(TpM::TangentSpace, X, V) + +Project the vector `V` from the embedding of the tangent space `TpM` (identified with ``T_X(T_p\mathcal M)``), +that is project the vector `V` onto the tangent space at `TpM.point`. +""" +project(::TangentSpace, ::Any, ::Any) + +function project!(TpM::TangentSpace, W, X, V) + return project!(TpM.manifold, W, TpM.point, V) +end + +function Random.rand!(TpM::TangentSpace, X; vector_at = nothing) + rand!(TpM.manifold, X; vector_at = TpM.point) + return X +end +function Random.rand!(rng::AbstractRNG, TpM::TangentSpace, X; vector_at = nothing) + rand!(rng, TpM.manifold, X; vector_at = TpM.point) + return X +end + +function representation_size(TpM::TangentSpace) + return representation_size(TpM.manifold) +end + +function Base.show(io::IO, ::MIME"text/plain", TpM::TangentSpace) + println(io, "Tangent space to the manifold $(base_manifold(TpM)) at point:") + pre = " " + sp = sprint(show, "text/plain", TpM.point; context = io, sizehint = 0) + sp = replace(sp, '\n' => "\n$(pre)") + return print(io, pre, sp) +end + + +function Base.show(io::IO, ::MIME"text/plain", cTpM::CotangentSpace) + println(io, "Cotangent space to the manifold $(base_manifold(cTpM)) at point:") + pre = " " + sp = sprint(show, "text/plain", cTpM.point; context = io, sizehint = 0) + sp = replace(sp, '\n' => "\n$(pre)") + return print(io, pre, sp) +end + +@doc raw""" + Y = Weingarten(TpM::TangentSpace, X, V, A) + Weingarten!(TpM::TangentSpace, Y, p, X, V) + +Compute the Weingarten map ``\mathcal W_X`` at `X` on the [`TangentSpace`](@ref) `TpM` with respect to the +tangent vector ``V \in T_p\mathcal M`` and the normal vector ``A \in N_p\mathcal M``. + +Since this a flat space by itself, the result is always the zero tangent vector. +""" +Weingarten(::TangentSpace, ::Any, ::Any, ::Any) + +Weingarten!(::TangentSpace, W, X, V, A) = fill!(W, 0) + +@doc raw""" + zero_vector(TpM::TangentSpace, X) + +Zero tangent vector at point `X` from the [`TangentSpace`](@ref) `TpM`, +that is the zero tangent vector at point `TpM.point`. +""" +zero_vector(::TangentSpace, ::Any...) + +function zero_vector!(M::TangentSpace, V, X) + return zero_vector!(M.manifold, V, M.point) +end diff --git a/src/VectorFiber.jl b/src/VectorFiber.jl new file mode 100644 index 00000000..3016e38d --- /dev/null +++ b/src/VectorFiber.jl @@ -0,0 +1,11 @@ + +""" + VectorSpaceFiber{𝔽,M,TSpaceType} = Fiber{𝔽,TSpaceType,M} + where {𝔽,M<:AbstractManifold{𝔽},TSpaceType<:VectorSpaceType} + +Alias for a [`Fiber`](@ref) when the fiber is a vector space. +""" +const VectorSpaceFiber{𝔽,M,TSpaceType} = + Fiber{𝔽,TSpaceType,M} where {𝔽,M<:AbstractManifold{𝔽},TSpaceType<:VectorSpaceType} + +LinearAlgebra.norm(M::VectorSpaceFiber, p, X) = norm(M.manifold, M.point, X) diff --git a/src/bases.jl b/src/bases.jl index 3d41edc3..68e4a7ed 100644 --- a/src/bases.jl +++ b/src/bases.jl @@ -23,17 +23,24 @@ Optionally: * [`representation_size`](@ref), * broadcasting for basic operations. """ -abstract type VectorSpaceType end +abstract type VectorSpaceType <: FiberType end +""" + struct TangentSpaceType <: VectorSpaceType end + +A type that indicates that a [`Fiber`](@ref) is a [`TangentSpace`](@ref). +""" struct TangentSpaceType <: VectorSpaceType end +""" + struct CotangentSpaceType <: VectorSpaceType end + +A type that indicates that a [`Fiber`](@ref) is a [`CotangentSpace`](@ref). +""" struct CotangentSpaceType <: VectorSpaceType end TCoTSpaceType = Union{TangentSpaceType,CotangentSpaceType} -const TangentSpace = TangentSpaceType() -const CotangentSpace = CotangentSpaceType() - """ AbstractBasis{𝔽,VST<:VectorSpaceType} @@ -65,14 +72,14 @@ for the vectors elements. struct DefaultBasis{𝔽,VST<:VectorSpaceType} <: AbstractBasis{𝔽,VST} vector_space::VST end -function DefaultBasis(𝔽::AbstractNumbers = ℝ, vs::VectorSpaceType = TangentSpace) +function DefaultBasis(𝔽::AbstractNumbers = ℝ, vs::VectorSpaceType = TangentSpaceType()) return DefaultBasis{𝔽,typeof(vs)}(vs) end -function DefaultBasis{𝔽}(vs::VectorSpaceType = TangentSpace) where {𝔽} +function DefaultBasis{𝔽}(vs::VectorSpaceType = TangentSpaceType()) where {𝔽} return DefaultBasis{𝔽,typeof(vs)}(vs) end function DefaultBasis{𝔽,TangentSpaceType}() where {𝔽} - return DefaultBasis{𝔽,TangentSpaceType}(TangentSpace) + return DefaultBasis{𝔽,TangentSpaceType}(TangentSpaceType()) end """ @@ -106,14 +113,17 @@ for the vectors elements. struct DefaultOrthogonalBasis{𝔽,VST<:VectorSpaceType} <: AbstractOrthogonalBasis{𝔽,VST} vector_space::VST end -function DefaultOrthogonalBasis(𝔽::AbstractNumbers = ℝ, vs::VectorSpaceType = TangentSpace) +function DefaultOrthogonalBasis( + 𝔽::AbstractNumbers = ℝ, + vs::VectorSpaceType = TangentSpaceType(), +) return DefaultOrthogonalBasis{𝔽,typeof(vs)}(vs) end -function DefaultOrthogonalBasis{𝔽}(vs::VectorSpaceType = TangentSpace) where {𝔽} +function DefaultOrthogonalBasis{𝔽}(vs::VectorSpaceType = TangentSpaceType()) where {𝔽} return DefaultOrthogonalBasis{𝔽,typeof(vs)}(vs) end function DefaultOrthogonalBasis{𝔽,TangentSpaceType}() where {𝔽} - return DefaultOrthogonalBasis{𝔽,TangentSpaceType}(TangentSpace) + return DefaultOrthogonalBasis{𝔽,TangentSpaceType}(TangentSpaceType()) end @@ -137,7 +147,7 @@ abstract type AbstractOrthonormalBasis{𝔽,VST<:VectorSpaceType} <: AbstractOrthogonalBasis{𝔽,VST} end """ - DefaultOrthonormalBasis(𝔽::AbstractNumbers = ℝ, vs::VectorSpaceType = TangentSpace) + DefaultOrthonormalBasis(𝔽::AbstractNumbers = ℝ, vs::VectorSpaceType = TangentSpaceType()) An arbitrary orthonormal basis of vector space of type `VST` on a manifold. This will usually be the fastest orthonormal basis available for a manifold. @@ -153,14 +163,17 @@ struct DefaultOrthonormalBasis{𝔽,VST<:VectorSpaceType} <: AbstractOrthonormal vector_space::VST end -function DefaultOrthonormalBasis(𝔽::AbstractNumbers = ℝ, vs::VectorSpaceType = TangentSpace) +function DefaultOrthonormalBasis( + 𝔽::AbstractNumbers = ℝ, + vs::VectorSpaceType = TangentSpaceType(), +) return DefaultOrthonormalBasis{𝔽,typeof(vs)}(vs) end -function DefaultOrthonormalBasis{𝔽}(vs::VectorSpaceType = TangentSpace) where {𝔽} +function DefaultOrthonormalBasis{𝔽}(vs::VectorSpaceType = TangentSpaceType()) where {𝔽} return DefaultOrthonormalBasis{𝔽,typeof(vs)}(vs) end function DefaultOrthonormalBasis{𝔽,TangentSpaceType}() where {𝔽} - return DefaultOrthonormalBasis{𝔽,TangentSpaceType}(TangentSpace) + return DefaultOrthonormalBasis{𝔽,TangentSpaceType}(TangentSpaceType()) end """ @@ -376,14 +389,14 @@ function _dual_basis( p, ::DefaultOrthonormalBasis{𝔽,TangentSpaceType}, ) where {𝔽} - return DefaultOrthonormalBasis{𝔽}(CotangentSpace) + return DefaultOrthonormalBasis{𝔽}(CotangentSpaceType()) end function _dual_basis( ::AbstractManifold, p, ::DefaultOrthonormalBasis{𝔽,CotangentSpaceType}, ) where {𝔽} - return DefaultOrthonormalBasis{𝔽}(TangentSpace) + return DefaultOrthonormalBasis{𝔽}(TangentSpaceType()) end function _euclidean_basis_vector(p::StridedArray, i) diff --git a/src/point_vector_fallbacks.jl b/src/point_vector_fallbacks.jl index 5b36dab5..797e67cb 100644 --- a/src/point_vector_fallbacks.jl +++ b/src/point_vector_fallbacks.jl @@ -331,6 +331,26 @@ macro default_manifold_fallbacks(TM, TP, TV, pfield::Symbol, vfield::Symbol) ManifoldsBase.retract_embedded!(M, q.$pfield, p.$pfield, X.$vfield, t, m) return q end + function ManifoldsBase.retract_sasaki( + M::$TM, + p::$TP, + X::$TV, + t::Number, + m::SasakiRetraction, + ) + return $TP(ManifoldsBase.retract_sasaki(M, p.$pfield, X.$vfield, t, m)) + end + function ManifoldsBase.retract_sasaki!( + M::$TM, + q::$TP, + p::$TP, + X::$TV, + t::Number, + m::SasakiRetraction, + ) + ManifoldsBase.retract_sasaki!(M, q.$pfield, p.$pfield, X.$vfield, t, m) + return q + end end, ) for f_postfix in [:polar, :project, :qr, :softmax] diff --git a/src/retractions.jl b/src/retractions.jl index 4892db13..6721511e 100644 --- a/src/retractions.jl +++ b/src/retractions.jl @@ -49,6 +49,7 @@ Retraction using the exponential map. """ struct ExponentialRetraction <: AbstractRetractionMethod end + @doc raw""" ODEExponentialRetraction{T<:AbstractRetractionMethod, B<:AbstractBasis} <: AbstractRetractionMethod @@ -173,6 +174,29 @@ function RetractionWithKeywords(m::T; kwargs...) where {T<:AbstractRetractionMet return RetractionWithKeywords{T,typeof(kwargs)}(m, kwargs) end +@doc raw""" + struct SasakiRetraction <: AbstractRetractionMethod end + +Exponential map on [`TangentBundle`](https://juliamanifolds.github.io/Manifolds.jl/stable/manifolds/vector_bundle.html#Manifolds.TangentBundle) computed via Euler integration as described +in [MuralidharanFletcher:2012](@cite). The system of equations for ``\gamma : ℝ \to T\mathcal M`` such that +``\gamma(1) = \exp_{p,X}(X_M, X_F)`` and ``\gamma(0)=(p, X)`` reads + +```math +\dot{\gamma}(t) = (\dot{p}(t), \dot{X}(t)) = (R(X(t), \dot{X}(t))\dot{p}(t), 0) +``` + +where ``R`` is the Riemann curvature tensor (see [`riemann_tensor`](@ref)). + +# Constructor + + SasakiRetraction(L::Int) + +In this constructor `L` is the number of integration steps. +""" +struct SasakiRetraction <: AbstractRetractionMethod + L::Int +end + """ SoftmaxRetraction <: AbstractRetractionMethod @@ -249,7 +273,7 @@ struct LogarithmicInverseRetraction <: AbstractInverseRetractionMethod end @doc raw""" PadeInverseRetraction{m} <: AbstractInverseRetractionMethod -An inverse retraction based on the Padé approximation of order $m$ for the retraction. +An inverse retraction based on the Padé approximation of order ``m`` for the retraction. !!! note "Technical Note" Though you would call e.g. [`inverse_retract`](@ref)`(M, p, q, PadeInverseRetraction(m))`, @@ -927,6 +951,9 @@ end function _retract!(M::AbstractManifold, q, p, X, t::Number, ::QRRetraction; kwargs...) return retract_qr!(M, q, p, X, t; kwargs...) end +function _retract!(M::AbstractManifold, q, p, X, t::Number, m::SasakiRetraction) + return retract_sasaki!(M, q, p, X, t, m) +end function _retract!(M::AbstractManifold, q, p, X, t::Number, ::SoftmaxRetraction; kwargs...) return retract_softmax!(M, q, p, X, t; kwargs...) end @@ -1047,6 +1074,15 @@ retract_softmax!(M::AbstractManifold, q, p, X, t::Number) function retract_softmax! end +""" + retract_sasaki!(M::AbstractManifold, q, p, X, t::Number, m::SasakiRetraction) + +Compute the in-place variant of the [`SasakiRetraction`](@ref) `m`. +""" +retract_sasaki!(M::AbstractManifold, q, p, X, t::Number, m::SasakiRetraction) + +function retract_sasaki! end + @doc raw""" retract(M::AbstractManifold, p, X, method::AbstractRetractionMethod=default_retraction_method(M, typeof(p))) retract(M::AbstractManifold, p, X, t::Number=1, method::AbstractRetractionMethod=default_retraction_method(M, typeof(p))) @@ -1112,6 +1148,9 @@ end function _retract(M::AbstractManifold, p, X, t::Number, ::QRRetraction; kwargs...) return retract_qr(M, p, X, t; kwargs...) end +function _retract(M::AbstractManifold, p, X, t::Number, m::SasakiRetraction) + return retract_sasaki(M, p, X, t, m) +end function _retract(M::AbstractManifold, p, X, t::Number, ::SoftmaxRetraction; kwargs...) return retract_softmax(M, p, X, t; kwargs...) end @@ -1238,5 +1277,16 @@ function retract_pade(M::AbstractManifold, p, X, t::Number, m::PadeRetraction; k return retract_pade!(M, q, p, X, t, m; kwargs...) end +""" + retract_sasaki(M::AbstractManifold, p, X, t::Number, m::SasakiRetraction) + +Compute the allocating variant of the [`SasakiRetraction`](@ref), +which by default allocates and calls `retract_sasaki!`. +""" +function retract_sasaki(M::AbstractManifold, p, X, t::Number, m::SasakiRetraction) + q = allocate_result(M, retract, p, X) + return retract_sasaki!(M, q, p, X, t, m) +end + Base.show(io::IO, ::CayleyRetraction) = print(io, "CayleyRetraction()") -Base.show(io::IO, ::PadeRetraction{m}) where {m} = print(io, "PadeRetraction($(m))") +Base.show(io::IO, ::PadeRetraction{m}) where {m} = print(io, "PadeRetraction($m)") diff --git a/src/vector_spaces.jl b/src/vector_spaces.jl index c00a498e..c1d845a2 100644 --- a/src/vector_spaces.jl +++ b/src/vector_spaces.jl @@ -53,10 +53,10 @@ const TFVector = FVector{TangentSpaceType} const CoTFVector = FVector{CotangentSpaceType} function TFVector(data, basis::AbstractBasis) - return TFVector{typeof(data),typeof(basis)}(TangentSpace, data, basis) + return TFVector{typeof(data),typeof(basis)}(TangentSpaceType(), data, basis) end function CoTFVector(data, basis::AbstractBasis) - return CoTFVector{typeof(data),typeof(basis)}(CotangentSpace, data, basis) + return CoTFVector{typeof(data),typeof(basis)}(CotangentSpaceType(), data, basis) end function Base.show(io::IO, fX::TFVector) @@ -125,9 +125,6 @@ function number_eltype( end number_eltype(v::FVector) = number_eltype(v.data) -Base.show(io::IO, ::TangentSpaceType) = print(io, "TangentSpace") -Base.show(io::IO, ::CotangentSpaceType) = print(io, "CotangentSpace") - """ vector_space_dimension(M::AbstractManifold, V::VectorSpaceType) diff --git a/src/vector_transport.jl b/src/vector_transport.jl index 42996b92..4d9dc197 100644 --- a/src/vector_transport.jl +++ b/src/vector_transport.jl @@ -36,7 +36,7 @@ by differentiation this retraction, is given by = D_Y\operatorname{retr}_p(Y)[X] = \frac{\mathrm{d}}{\mathrm{d}t}\operatorname{retr}_p(Y+tX)\Bigr|_{t=0}. ``` -see [^AbsilMahonySepulchre2008], Section 8.1.2 for more details. +see [AbsilMahonySepulchre:2008](@cite), Section 8.1.2 for more details. This can be phrased similarly as a [`vector_transport_to`](@ref) by introducing ``q=\operatorname{retr}_pX`` and defining @@ -92,9 +92,9 @@ The [`PoleLadderTransport`](@ref) posesses two advantages compared to * it is cheaper to evaluate, if you want to transport several vectors, since the mid point $c$ then stays unchanged. * while both methods are exact if the curvature is zero, pole ladder is even exact in - symmetric Riemannian manifolds[^Pennec2018] + symmetric Riemannian manifolds [Pennec:2018](@cite) -The pole ladder was was proposed in [^LorenziPennec2014]. Its name stems from the fact that +The pole ladder was was proposed in [LorenziPennec:2013](@cite). Its name stems from the fact that it resembles a pole ladder when applied to a sequence of points usccessively. # Constructor @@ -105,20 +105,9 @@ PoleLadderTransport( ) ```` Construct the classical pole ladder that employs exp and log, i.e. as proposed -in[^LorenziPennec2014]. For an even cheaper transport the inner operations can be +in[LorenziPennec:2013](@cite). For an even cheaper transport the inner operations can be changed to an [`AbstractRetractionMethod`](@ref) `retraction` and an [`AbstractInverseRetractionMethod`](@ref) `inverse_retraction`, respectively. - -[^LorenziPennec2014]: - > Lorenzi, M. and Pennec, X: Efficient parallel transport of deformations in time - > series of images: From Schild’s to pole ladder. - > Journal of Mathematical Imaging and Vision (2014), 50(1), pp. 5–17 - > doi [10.1007/s10851-013-0470-3](https://doi.org/10.1007/s10851-013-0470-3), - > hal: [hal-00870489](https://hal.inria.fr/hal-00870489) -[^Pennec2018]: - > Pennec, X: Parallel Transport with Pole Ladder: a Third Order Scheme in Affine - > Connection Spaces which is Exact in Affine Symmetric Spaces. - > arXiv: [1805.11436](https://arxiv.org/abs/1805.11436) """ struct PoleLadderTransport{ RT<:AbstractRetractionMethod, @@ -141,7 +130,7 @@ end ScaledVectorTransport{T} <: AbstractVectorTransportMethod Introduce a scaled variant of any [`AbstractVectorTransportMethod`](@ref) `T`, -as introduced in [^SatoIwai2013] for some ``X∈ T_p\mathcal M`` as +as introduced in [SatoIwai:2013](@cite) for some ``X∈ T_p\mathcal M`` as ```math \mathcal T^{\mathrm{S}}(X) = \frac{\lVert X\rVert_p}{\lVert \mathcal T(X)\rVert_q}\mathcal T(X). @@ -154,12 +143,6 @@ retraction). Therefore a default implementation is only provided for the [`vecto # Constructor ScaledVectorTransport(m::AbstractVectorTransportMethod) - -[^SatoIwai2013]: - > Sato, H., Iwai, T.: _A new, globally convergent Riemannian conjugate gradient method_, - > Optimization, 2013, Volume 64(4), pp. 1011–1031. - > doi: [10.1080/02331934.2013.836650](https://doi.org/10.1080/02331934.2013.836650), - > arXiv: [1302.0125](https://arxiv.org/abs/1302.0125). """ struct ScaledVectorTransport{T<:AbstractVectorTransportMethod} <: AbstractVectorTransportMethod @@ -187,7 +170,7 @@ This method employs the internal function [`schilds_ladder`](@ref)`(M, p, d, q)` leaving the manifold. The name stems from the image of this paralleltogram in a repeated application yielding the -image of a ladder. The approximation was proposed in [^EhlersPiraniSchild1972]. +image of a ladder. The approximation was proposed in [EhlersPiraniSchild:1972](@cite). # Constructor ````julia @@ -197,15 +180,9 @@ SchildsLadderTransport( ) ```` Construct the classical Schilds ladder that employs exp and log, i.e. as proposed -in[^EhlersPiraniSchild1972]. For an even cheaper transport these inner operations can be +in [EhlersPiraniSchild:1972](@cite). For an even cheaper transport these inner operations can be changed to an [`AbstractRetractionMethod`](@ref) `retraction` and an [`AbstractInverseRetractionMethod`](@ref) `inverse_retraction`, respectively. - -[^EhlersPiraniSchild1972]: - > Ehlers, J., Pirani, F.A.E., Schild, A.: The geometry of free fall and light - > propagation. In: O’Raifeartaigh, L. (ed.) General Relativity: Papers in Honour of - > J. L. Synge, pp. 63–84. Clarendon Press, Oxford (1972). - > reprint doi: [10.1007/s10714-012-1353-4](https://doi.org/10.1007/s10714-012-1353-4) """ struct SchildsLadderTransport{ RT<:AbstractRetractionMethod, @@ -875,7 +852,7 @@ end Given an [`AbstractManifold`](@ref) ``\mathcal M`` the vector transport is a generalization of the [`parallel_transport_direction`](@ref) that identifies vectors from different tangent spaces. -More precisely using [^AbsilMahonySepulchre2008], Def. 8.1.1, a vector transport +More precisely using [AbsilMahonySepulchre:2008](@cite), Def. 8.1.1, a vector transport ``T_{p,d}: T_p\mathcal M \to T_q\mathcal M``, ``p∈ \mathcal M``, ``Y∈ T_p\mathcal M`` is a smooth mapping associated to a retraction ``\operatorname{retr}_p(Y) = q`` such that @@ -903,13 +880,6 @@ Instead of spcifying a start direction `d` one can equivalently also specify a t ``T_q\mathcal M``, see [`vector_transport_to`](@ref). By default [`vector_transport_direction`](@ref) falls back to using [`vector_transport_to`](@ref), using the [`default_retraction_method`](@ref) on `M`. - -[^AbsilMahonySepulchre2008]: - > Absil, P.-A., Mahony, R. and Sepulchre R., - > _Optimization Algorithms on Matrix Manifolds_ - > Princeton University Press, 2008, - > doi: [10.1515/9781400830244](https://doi.org/10.1515/9781400830244) - > [open access](http://press.princeton.edu/chapters/absil/) """ function vector_transport_direction( M::AbstractManifold, diff --git a/test/bases.jl b/test/bases.jl index ad5e64ec..efbaa935 100644 --- a/test/bases.jl +++ b/test/bases.jl @@ -1,7 +1,7 @@ using LinearAlgebra using ManifoldsBase using ManifoldsBase: DefaultManifold, ℝ, ℂ, RealNumbers, ComplexNumbers -using ManifoldsBase: CotangentSpace, CotangentSpaceType, TangentSpace, TangentSpaceType +using ManifoldsBase: CotangentSpaceType, TangentSpaceType using ManifoldsBase: FVector using Test import Base: +, -, *, copyto!, isapprox @@ -267,11 +267,12 @@ DiagonalizingBasisProxy() = DiagonalizingOrthonormalBasis([1.0, 0.0, 0.0]) @test DefaultOrthonormalBasis{ℂ,TangentSpaceType}() === DefaultOrthonormalBasis(ℂ) - @test DefaultBasis{ℂ}(CotangentSpace) === DefaultBasis(ℂ, CotangentSpace) - @test DefaultOrthogonalBasis{ℂ}(CotangentSpace) === - DefaultOrthogonalBasis(ℂ, CotangentSpace) - @test DefaultOrthonormalBasis{ℂ}(CotangentSpace) === - DefaultOrthonormalBasis(ℂ, CotangentSpace) + @test DefaultBasis{ℂ}(CotangentSpaceType()) === + DefaultBasis(ℂ, CotangentSpaceType()) + @test DefaultOrthogonalBasis{ℂ}(CotangentSpaceType()) === + DefaultOrthogonalBasis(ℂ, CotangentSpaceType()) + @test DefaultOrthonormalBasis{ℂ}(CotangentSpaceType()) === + DefaultOrthonormalBasis(ℂ, CotangentSpaceType()) end _pts = [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]] @@ -403,7 +404,7 @@ DiagonalizingBasisProxy() = DiagonalizingOrthonormalBasis([1.0, 0.0, 0.0]) M = DefaultManifold(2, 3) x = collect(reshape(1.0:6.0, (2, 3))) pb = get_basis(M, x, DefaultOrthonormalBasis()) - B2 = DefaultOrthonormalBasis(ManifoldsBase.ℝ, ManifoldsBase.CotangentSpace) + B2 = DefaultOrthonormalBasis(ManifoldsBase.ℝ, ManifoldsBase.CotangentSpaceType()) pb2 = get_basis(M, x, B2) test_basis_string = """ @@ -496,24 +497,24 @@ DiagonalizingBasisProxy() = DiagonalizingOrthonormalBasis([1.0, 0.0, 0.0]) end @testset "Bases of cotangent spaces" begin - b1 = DefaultOrthonormalBasis(ℝ, CotangentSpace) - @test b1.vector_space == CotangentSpace + b1 = DefaultOrthonormalBasis(ℝ, CotangentSpaceType()) + @test b1.vector_space == CotangentSpaceType() - b2 = DefaultOrthogonalBasis(ℝ, CotangentSpace) - @test b2.vector_space == CotangentSpace + b2 = DefaultOrthogonalBasis(ℝ, CotangentSpaceType()) + @test b2.vector_space == CotangentSpaceType() - b3 = DefaultBasis(ℝ, CotangentSpace) - @test b3.vector_space == CotangentSpace + b3 = DefaultBasis(ℝ, CotangentSpaceType()) + @test b3.vector_space == CotangentSpaceType() M = DefaultManifold(2; field = ℂ) p = [1.0, 2.0im] b1_d = ManifoldsBase.dual_basis(M, p, b1) @test b1_d isa DefaultOrthonormalBasis - @test b1_d.vector_space == TangentSpace + @test b1_d.vector_space == TangentSpaceType() b1_d_d = ManifoldsBase.dual_basis(M, p, b1_d) @test b1_d_d isa DefaultOrthonormalBasis - @test b1_d_d.vector_space == CotangentSpace + @test b1_d_d.vector_space == CotangentSpaceType() end @testset "Complex Basis - Mutating cases" begin @@ -544,25 +545,25 @@ DiagonalizingBasisProxy() = DiagonalizingOrthonormalBasis([1.0, 0.0, 0.0]) end @testset "FVector" begin - @test sprint(show, TangentSpace) == "TangentSpace" - @test sprint(show, CotangentSpace) == "CotangentSpace" + @test sprint(show, TangentSpaceType()) == "TangentSpaceType()" + @test sprint(show, CotangentSpaceType()) == "CotangentSpaceType()" tvs = ([1.0, 0.0, 0.0], [0.0, 1.0, 0.0]) fv_tvs = map(v -> TFVector(v, DefaultOrthonormalBasis()), tvs) fv1 = fv_tvs[1] tv1s = allocate(fv_tvs[1]) @test isa(tv1s, FVector) - @test tv1s.type == TangentSpace + @test tv1s.type == TangentSpaceType() @test size(tv1s.data) == size(tvs[1]) @test number_eltype(tv1s) == number_eltype(tvs[1]) @test number_eltype(tv1s) == number_eltype(typeof(tv1s)) @test isa(fv1 + fv1, FVector) - @test (fv1 + fv1).type == TangentSpace + @test (fv1 + fv1).type == TangentSpaceType() @test isa(fv1 - fv1, FVector) - @test (fv1 - fv1).type == TangentSpace + @test (fv1 - fv1).type == TangentSpaceType() @test isa(-fv1, FVector) - @test (-fv1).type == TangentSpace + @test (-fv1).type == TangentSpaceType() @test isa(2 * fv1, FVector) - @test (2 * fv1).type == TangentSpace + @test (2 * fv1).type == TangentSpaceType() tv1s_32 = allocate(fv_tvs[1], Float32) @test isa(tv1s, FVector) @test eltype(tv1s_32.data) === Float32 @@ -571,7 +572,7 @@ DiagonalizingBasisProxy() = DiagonalizingOrthonormalBasis([1.0, 0.0, 0.0]) @test sprint(show, fv1) == "TFVector([1.0, 0.0, 0.0], $(fv1.basis))" - cofv1 = CoTFVector(tvs[1], DefaultOrthonormalBasis(ℝ, CotangentSpace)) + cofv1 = CoTFVector(tvs[1], DefaultOrthonormalBasis(ℝ, CotangentSpaceType())) @test cofv1 isa CoTFVector @test sprint(show, cofv1) == "CoTFVector([1.0, 0.0, 0.0], $(fv1.basis))" end @@ -579,10 +580,10 @@ DiagonalizingBasisProxy() = DiagonalizingOrthonormalBasis([1.0, 0.0, 0.0]) @testset "vector_space_dimension" begin M = ManifoldsBase.DefaultManifold(3) MC = ManifoldsBase.DefaultManifold(3; field = ℂ) - @test ManifoldsBase.vector_space_dimension(M, TangentSpace) == 3 - @test ManifoldsBase.vector_space_dimension(M, CotangentSpace) == 3 - @test ManifoldsBase.vector_space_dimension(MC, TangentSpace) == 6 - @test ManifoldsBase.vector_space_dimension(MC, CotangentSpace) == 6 + @test ManifoldsBase.vector_space_dimension(M, TangentSpaceType()) == 3 + @test ManifoldsBase.vector_space_dimension(M, CotangentSpaceType()) == 3 + @test ManifoldsBase.vector_space_dimension(MC, TangentSpaceType()) == 6 + @test ManifoldsBase.vector_space_dimension(MC, CotangentSpaceType()) == 6 end @testset "requires_caching" begin diff --git a/test/default_manifold.jl b/test/default_manifold.jl index 73e34b3f..4c8c1e22 100644 --- a/test/default_manifold.jl +++ b/test/default_manifold.jl @@ -181,6 +181,16 @@ function ManifoldsBase.retract_exp_ode!( return (q .= p .+ t .* X) end ManifoldsBase.retract_pade!(::DefaultManifold, q, p, X, t::Number, i) = (q .= p .+ t .* X) +function ManifoldsBase.retract_sasaki!( + ::DefaultManifold, + q, + p, + X, + t::Number, + ::SasakiRetraction, +) + return (q .= p .+ t .* X) +end ManifoldsBase.retract_softmax!(::DefaultManifold, q, p, X, t::Number) = (q .= p .+ t .* X) ManifoldsBase.get_embedding(M::DefaultManifold) = M # dummy embedding ManifoldsBase.inverse_retract_polar!(::DefaultManifold, Y, p, q) = (Y .= q .- p) @@ -723,6 +733,7 @@ Base.size(x::MatrixVectorTransport) = (size(x.m, 2),) ODEExponentialRetraction(PolarRetraction(), DefaultBasis()), PadeRetraction(2), EmbeddedRetraction(ExponentialRetraction()), + SasakiRetraction(5), ] @test retract(M, q, Y, retr) == DefaultPoint(q.value + Y.value) @test retract(M, q, Y, 0.5, retr) == DefaultPoint(q.value + 0.5 * Y.value) diff --git a/test/fibers.jl b/test/fibers.jl new file mode 100644 index 00000000..e64bf1f9 --- /dev/null +++ b/test/fibers.jl @@ -0,0 +1,101 @@ +using RecursiveArrayTools, ManifoldsBase, Test +using Random +using ManifoldsBase: DefaultManifold, VectorSpaceType, ℝ, Fiber +struct TestVectorSpaceType <: VectorSpaceType end + +struct TestFiberType <: ManifoldsBase.FiberType end + +function ManifoldsBase.fiber_dimension(M::AbstractManifold, ::TestFiberType) + return 2 * manifold_dimension(M) +end + +function ManifoldsBase.vector_space_dimension(M::AbstractManifold, ::TestVectorSpaceType) + return 2 * manifold_dimension(M) +end + +include("test_manifolds.jl") + + +@testset "vector space fibers" begin + M = DefaultManifold(3) + + p = [1.0, 0.0, 0.0] + + TpM = TangentSpace(M, p) + @test is_flat(TpM) + + @test ManifoldsBase.fiber_dimension(M, CotangentSpaceType()) == 3 + + @test ManifoldsBase.fiber_dimension(M, ManifoldsBase.TangentSpaceType()) == 3 + + @testset "spaces at point" begin + p = [1.0, 0.0, 0.0] + q = [0.0, 2.0, 0.0] + t_p = TangentSpace(M, p) + t_ps = sprint(show, "text/plain", t_p) + sp = sprint(show, "text/plain", p) + sp = replace(sp, '\n' => "\n ") + t_ps_test = "Tangent space to the manifold $(M) at point:\n $(sp)" + @test t_ps == t_ps_test + ct_p = CotangentSpace(M, p) + ct_ps = sprint(show, "text/plain", ct_p) + ct_ps_test = "Cotangent space to the manifold $(M) at point:\n $(sp)" + @test ct_ps == ct_ps_test + @test base_manifold(t_p) == M + @test manifold_dimension(t_p) == 3 + @test t_p.manifold == M + @test t_p.fiber_type == TangentSpaceType() + @test t_p.point == p + @test injectivity_radius(t_p) == Inf + @test representation_size(t_p) == representation_size(M) + X = [0.0, 0.0, 1.0] + Y = [0.0, 2.0, -1.0] + @test embed(t_p, X) == X + @test embed(t_p, X, X) == X + @test distance(t_p, p, q) ≈ sqrt(5) + @test isapprox(t_p, exp(t_p, X, X), 2 * X) + @test isapprox(t_p, log(t_p, X, Y), [0.0, 2.0, -2.0]) + @test isapprox(t_p, X, log(t_p, X, Y), [0.0, 2.0, -2.0]) + @test inner(t_p, X, X, X) ≈ 1.0 + @test norm(t_p, X, X) ≈ 1.0 + @test parallel_transport_to(t_p, X, Y, X) ≈ Y + @test vector_transport_to(t_p, X, Y, X) ≈ Y + @test vector_transport_to(t_p, X, Y, X, ProjectionTransport()) ≈ Y + Z = similar(X) + @test vector_transport_to!(t_p, Z, X, Y, X, ProjectionTransport()) ≈ Y + @test project(t_p, X, Y) ≈ Y + @test project(t_p, Y) ≈ Y + @test rand(t_p) isa Vector{Float64} + @test rand(t_p; vector_at = X) isa Vector{Float64} + @test rand(Random.default_rng(), t_p) isa Vector{Float64} + @test rand(Random.default_rng(), t_p; vector_at = X) isa Vector{Float64} + # generic vector space at + X_p = Fiber(M, p, TestVectorSpaceType()) + X_ps = sprint(show, "text/plain", X_p) + X_ps_test = "VectorSpaceFiber{ℝ, DefaultManifold{ℝ, Tuple{Int64}}, TestVectorSpaceType, Vector{Float64}}\nFiber:\n TestVectorSpaceType()DefaultManifold(3; field = ℝ)\nBase point:\n $(sp)" + @test X_ps == X_ps_test + + for basis in + [DefaultOrthonormalBasis(), get_basis(t_p, p, DefaultOrthonormalBasis())] + @test length(get_vectors(t_p, p, get_basis(t_p, p, basis))) == 3 + X1c = get_coordinates(t_p, p, X, basis) + @test isapprox(X1c, [0.0, 0.0, 1.0]) + Y1c = similar(X1c) + get_coordinates!(t_p, Y1c, p, X, basis) + @test isapprox(X1c, Y1c) + @test isapprox(get_vector(t_p, p, X1c, basis), X) + Z1 = similar(X) + get_vector!(t_p, Z1, p, X1c, basis) + @test isapprox(Z1, X) + end + end + + @testset "Weingarten Map" begin + p0 = [1.0, 0.0, 0.0] + M = TangentSpace(M, p0) + p = [0.0, 1.0, 1.0] + X = [0.0, 1.0, 0.0] + V = [1.0, 0.0, 0.0] + @test Weingarten(M, p, X, V) == zero_vector(M, p) + end +end diff --git a/test/power.jl b/test/power.jl index b5494922..87192833 100644 --- a/test/power.jl +++ b/test/power.jl @@ -1,10 +1,15 @@ using Test using ManifoldsBase -using ManifoldsBase: AbstractNumbers, ℝ, ℂ, NestedReplacingPowerRepresentation +using ManifoldsBase: + AbstractNumbers, ℝ, ℂ, NestedReplacingPowerRepresentation, VectorSpaceType using StaticArrays using LinearAlgebra using Random +include("test_manifolds.jl") + +struct TestVectorSpaceType <: VectorSpaceType end + power_array_wrapper(::Type{NestedPowerRepresentation}, ::Int) = identity power_array_wrapper(::Type{NestedReplacingPowerRepresentation}, i::Int) = SVector{i} @@ -293,4 +298,60 @@ struct TestArrayRepresentation <: AbstractPowerRepresentation end rand(MersenneTwister(123), N; vector_at = p) end end + + @testset "metric conversion" begin + M = TestSPD(3) + N = PowerManifold(M, NestedPowerRepresentation(), 2) + e = EuclideanMetric() + p = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1] + q = [2.0 0.0 0.0; 0.0 2.0 0.0; 0.0 0.0 1] + P = [p, q] + X = [log(M, p, q), log(M, q, p)] + Y = change_metric(N, e, P, X) + Yc = [change_metric(M, e, p, log(M, p, q)), change_metric(M, e, q, log(M, q, p))] + @test norm(N, P, Y .- Yc) ≈ 0 + Z = change_representer(N, e, P, X) + Zc = [ + change_representer(M, e, p, log(M, p, q)), + change_representer(M, e, q, log(M, q, p)), + ] + @test norm(N, P, Z .- Zc) ≈ 0 + end + + @testset "Other stuff" begin + M1 = TestSphere(2) + @testset "Weingarten" begin + Mpr = PowerManifold(M1, NestedPowerRepresentation(), 2) + p = [1.0, 0.0, 0.0] + X = [0.0, 0.2, 0.0] + V = [0.1, 0.0, 0.0] #orthogonal to TpM -> parallel to p + @test isapprox( + Mpr, + Weingarten(Mpr, [p, p], [X, X], [V, V]), + [-0.1 * X, -0.1 * X], + ) + end + end + + @testset "Type size" begin + M = ManifoldsBase.DefaultManifold(3) + N = PowerManifold(M, NestedPowerRepresentation(), 2; parameter = :type) + p = [[1, 2, 3], [3, 4, 5]] + @test zero_vector(N, p) == 0 .* p + @test PowerManifold(N, 3) isa PowerManifold{ + ℝ, + <:DefaultManifold, + ManifoldsBase.TypeParameter{Tuple{2,3}}, + NestedPowerRepresentation, + } + @test PowerManifold(N, 3; parameter = :field) isa + PowerManifold{ℝ,<:DefaultManifold,Tuple{Int,Int},NestedPowerRepresentation} + @test repr(N) == + "PowerManifold(DefaultManifold(3; field = ℝ), NestedPowerRepresentation(), 2; parameter=:type)" + + M = ManifoldsBase.DefaultManifold(3) + N = PowerManifold(M, NestedPowerRepresentation(), 2, 3; parameter = :type) + p = rand(N) + @test zero_vector(N, p) == 0 .* p + end end diff --git a/test/product_manifold.jl b/test/product_manifold.jl new file mode 100644 index 00000000..3f124c5c --- /dev/null +++ b/test/product_manifold.jl @@ -0,0 +1,560 @@ +using Test +using ManifoldsBase +using ManifoldsBase: submanifold_component, submanifold_components +using ManifoldsBase: + AbstractNumbers, ℝ, ℂ, NestedReplacingPowerRepresentation, ProductBasisData +using LinearAlgebra +using Random +using RecursiveArrayTools + +include("test_sphere.jl") + +@testset "Product manifold" begin + M1 = TestSphere(2) + M2 = ManifoldsBase.DefaultManifold(2, 2) + @test (@inferred ProductManifold(M1, M2)) isa ProductManifold + + M = ProductManifold(M1, M2) + + @test_throws MethodError ProductManifold() + + p1 = ArrayPartition([1, 0.0, 0.0], [4.0 5.0; 6.0 7.0]) + p2 = ArrayPartition([0.0, 1.0, 0.0], [4.0 8.0; 3.0 7.5]) + X1 = ArrayPartition([0.0, 1.0, 0.2], [4.0 0.0; 2.0 7.0]) + + @test !is_flat(M) + @test M[1] == M1 + @test M[2] == M2 + @test injectivity_radius(M) ≈ π + @test injectivity_radius( + M, + ProductRetraction(ExponentialRetraction(), ExponentialRetraction()), + ) ≈ π + @test injectivity_radius(M, ExponentialRetraction()) ≈ π + @test injectivity_radius( + M, + p1, + ProductRetraction(ExponentialRetraction(), ExponentialRetraction()), + ) ≈ π + @test injectivity_radius(M, p1) ≈ π + @test injectivity_radius(M, p1, ExponentialRetraction()) ≈ π + + @test ManifoldsBase.number_of_components(M) == 2 + @test ManifoldsBase.submanifold(M, 1) === M1 + @test ManifoldsBase.submanifold(M, [2, 1]) === ProductManifold(M2, M1) + @test ManifoldsBase.check_point(M, p1) === nothing + @test ManifoldsBase.check_vector(M, p1, X1) === nothing + # test that arrays are not points + @test_throws DomainError is_point(M, [1, 2], true) + @test ManifoldsBase.check_point(M, [1, 2]) isa DomainError + @test_throws DomainError is_vector(M, 1, [1, 2], true; check_base_point = false) + @test ManifoldsBase.check_vector(M, 1, [1, 2]; check_base_point = false) isa DomainError + #default fallbacks for check_size, Product not working with Arrays + @test ManifoldsBase.check_size(M, zeros(2)) isa DomainError + @test ManifoldsBase.check_size(M, zeros(2), zeros(3)) isa DomainError + + retraction_methods = [ + ManifoldsBase.ProductRetraction( + ManifoldsBase.ExponentialRetraction(), + ManifoldsBase.ExponentialRetraction(), + ), + ] + inverse_retraction_methods = [ + ManifoldsBase.InverseProductRetraction( + ManifoldsBase.LogarithmicInverseRetraction(), + ManifoldsBase.LogarithmicInverseRetraction(), + ), + ] + + @test ManifoldsBase.default_retraction_method(M) === retraction_methods[1] + @test ManifoldsBase.default_inverse_retraction_method(M) === + inverse_retraction_methods[1] + @test ManifoldsBase.default_inverse_retraction_method(M, typeof(X1)) === + inverse_retraction_methods[1] + + @testset "get_component, set_component!, getindex and setindex!" begin + @test get_component(M, p1, 1) == p1.x[1] + @test get_component(M, p1, Val(1)) == p1.x[1] + @test p1[M, 1] == p1.x[1] + @test p1[M, Val(1)] == p1.x[1] + @test p1[M, 1] isa Vector + @test p1[M, Val(1)] isa Vector + p2c = [5 6; 4 0] + set_component!(M, p1, p2c, 2) + @test get_component(M, p1, 2) == p2c + p1[M, 2] = 2 * p2c + @test p1[M, 2] == 2 * p2c + p3 = [11.0 15.0; 3 3] + set_component!(M, p1, p3, Val(2)) + @test get_component(M, p1, Val(2)) == p3 + p1[M, Val(2)] = 2 * p3 + @test p1[M, Val(2)] == 2 * p3 + + p1c = copy(p1) + p1c.x[1][1] = -123.0 + @test p1c.x[1][1] == -123.0 + @test p1.x[1][1] == 1.0 + copyto!(p1c, p1) + @test p1c.x[1][1] == 1.0 + + p1c.x[1][1] = -123.0 + copyto!(p1, p1c) + @test p1.x[1][1] == -123.0 + end + + @testset "some ArrayPartition functions" begin + q = allocate(p1) + @test q.x[1] isa Vector + p = ArrayPartition([[0.0, 1.0, 0.0]], [0.0, 0.0]) + q = allocate(p, Int) + @test q.x[1] isa Vector{Vector{Int}} + + q = allocate(p1) + copyto!(M, q, p1) + @test isapprox(q, p1) + @test ManifoldsBase.allocate_result(M, zero_vector) isa ArrayPartition + end + + @testset "allocate on PowerManifold of ProductManifold" begin + q = allocate([p1]) + @test q[1] isa ArrayPartition + @test q[1].x[1] isa Vector + end + + p1 = ArrayPartition([1.0, 0.0, 0.0], [4.0 5.0; 6.0 7.0]) + p2 = ArrayPartition([0.0, 1.0, 0.0], [4.0 8.0; 3.0 7.5]) + X1 = ArrayPartition([0.0, 1.0, 0.2], [4.0 0.0; 2.0 7.0]) + X2 = ArrayPartition([0.0, -1.0, 0.4], [3.0 1.0; -2.0 2.0]) + + @testset "Basic operations" begin + @test manifold_dimension(M) == 6 + @test representation_size(M) == (7,) + @test distance(M, p1, p2) ≈ 4.551637188998299 + qr = similar(p1) + exp!(M, qr, p1, X1) + @test exp(M, p1, X1) ≈ ArrayPartition( + [0.5235330372543839, 0.8354600062374664, 0.1670920012474933], + [8.0 5.0; 8.0 14.0], + ) + @test exp(M, p1, X1) ≈ qr + @test exp(M, p1, X1, 2.0) ≈ exp(M, p1, 2 * X1) + exp!(M, qr, p1, X1, 2.0) + @test qr ≈ exp(M, p1, 2 * X1) + @test qr ≈ retract( + M, + p1, + 2 * X1, + ProductRetraction(ExponentialRetraction(), ExponentialRetraction()), + ) + # single retraction gets “broadcasted” + @test qr ≈ retract(M, p1, 2 * X1, ExponentialRetraction()) + qr2 = similar(p1) + retract!( + M, + qr2, + p1, + 2 * X1, + ProductRetraction(ExponentialRetraction(), ExponentialRetraction()), + ) + @test qr2 ≈ qr + qr3 = similar(p1) + retract!(M, qr3, p1, 2 * X1, ExponentialRetraction()) + @test qr3 ≈ qr + Xr = similar(X1) + log!(M, Xr, p1, p2) + @test log(M, p1, p2) ≈ + ArrayPartition([0.0, 1.5707963267948966, 0.0], [0.0 3.0; -3.0 0.5]) + @test log(M, p1, p2) ≈ Xr + @test inverse_retract( + M, + p1, + p2, + InverseProductRetraction( + LogarithmicInverseRetraction(), + LogarithmicInverseRetraction(), + ), + ) ≈ Xr + # single inverse retraction gets “broadcasted” + @test inverse_retract(M, p1, p2, LogarithmicInverseRetraction()) ≈ Xr + + Zr = similar(X1) + inverse_retract!( + M, + Zr, + p1, + p2, + InverseProductRetraction( + LogarithmicInverseRetraction(), + LogarithmicInverseRetraction(), + ), + ) + @test Zr ≈ Xr + + @test inner(M, p1, X1, X2) ≈ 21.08 + @test norm(M, p1, X1) ≈ sqrt(70.04) + @test project(M, p1) ≈ p1 + q1 = similar(p1) + project!(M, q1, p1) + @test q1 ≈ p1 + @test project(M, p1, X1) ≈ X1 + Y1 = similar(X1) + project!(M, Y1, p1, X1) + @test X1 ≈ Y1 + @test zero_vector(M, p1) == zero(X1) + + @test rand(M) isa ArrayPartition + @test rand(M; vector_at = p1) isa ArrayPartition + @test rand(Random.default_rng(), M) isa ArrayPartition + @test rand(Random.default_rng(), M; vector_at = p1) isa ArrayPartition + @test is_point(M, rand!(M, q1)) + @test is_vector(M, p1, rand!(M, Y1; vector_at = p1)) + end + + @testset "Broadcasting" begin + br_result = p1 .+ 2.0 .* p2 + @test br_result isa ArrayPartition + @test br_result.x[1] ≈ [1.0, 2.0, 0.0] + @test br_result.x[2] ≈ [12.0 21.0; 12.0 22.0] + + br_result .= 2.0 .* p1 .+ p2 + @test br_result.x[1] ≈ [2.0, 1.0, 0.0] + @test br_result.x[2] ≈ [12.0 18.0; 15.0 21.5] + + br_result .= p1 + @test br_result.x[1] ≈ p1.x[1] + @test br_result.x[2] ≈ p1.x[2] + + @test axes(p1) == (Base.OneTo(7),) + + @test ManifoldsBase._get_vector_cache_broadcast(p1) === Val(false) + + # errors + p3 = ArrayPartition([3.0, 4.0, 5.0], [2.0, 5.0], [3.0, 2.0]) + @test_throws DimensionMismatch p1 .+ p3 + @test_throws DimensionMismatch p1 .= p3 + end + + @testset "CompositeManifoldError" begin + Mpr = ProductManifold(M1, M1) + let p1 = [1.0, 0.0, 0.0], + p2 = [0.0, 1.0, 0.0], + X1 = [0.0, 1.0, 0.2], + X2 = [1.0, 0.0, 0.2] + + p = ArrayPartition(p1, p2) + X = ArrayPartition(X1, X2) + pf = ArrayPartition(p1, X1) + Xf = ArrayPartition(X1, p2) + @test is_point(Mpr, p, true) + @test_throws CompositeManifoldError is_point(Mpr, X, true) + @test_throws ComponentManifoldError is_vector(Mpr, pf, X, true) + @test_throws ComponentManifoldError is_vector(Mpr, p, Xf, true) + end + end + + @testset "Weingarten" begin + Mpr = ProductManifold(M1, M1) + p = [1.0, 0.0, 0.0] + X = [0.0, 0.2, 0.0] + V = [0.1, 0.0, 0.0] #orthogonal to TpM -> parallel to p + @test isapprox( + M, + Weingarten( + Mpr, + ArrayPartition(p, p), + ArrayPartition(X, X), + ArrayPartition(V, V), + ), + ArrayPartition(-0.1 * X, -0.1 * X), + ) + end + + @testset "arithmetic" begin + @test isapprox(M, p1 + p2, ArrayPartition([1.0, 1.0, 0.0], [8.0 13.0; 9.0 14.5])) + @test isapprox(M, p1 - p2, ArrayPartition([1.0, -1.0, 0.0], [0.0 -3.0; 3.0 -0.5])) + @test isapprox(M, -p1, ArrayPartition([-1.0, -0.0, -0.0], [-4.0 -5.0; -6.0 -7.0])) + @test isapprox(M, p1 * 2, ArrayPartition([2.0, 0.0, 0.0], [8.0 10.0; 12.0 14.0])) + @test isapprox(M, 2 * p1, ArrayPartition([2.0, 0.0, 0.0], [8.0 10.0; 12.0 14.0])) + @test isapprox(M, p1 / 2, ArrayPartition([0.5, 0.0, 0.0], [2.0 2.5; 3.0 3.5])) + end + + @testset "Show methods" begin + M2p = ProductManifold(M1, M1, M2, M2) + @test sprint(show, M2p) == "ProductManifold($(M1), $(M1), $(M2), $(M2))" + withenv("LINES" => 10, "COLUMNS" => 100) do + @test sprint(show, "text/plain", ProductManifold(M1)) == + "ProductManifold with 1 submanifold:\n $(M1)" + @test sprint(show, "text/plain", M2p) == + "ProductManifold with 4 submanifolds:\n $(M1)\n $(M1)\n $(M2)\n $(M2)" + return nothing + end + withenv("LINES" => 7, "COLUMNS" => 100) do + @test sprint(show, "text/plain", M2p) == + "ProductManifold with 4 submanifolds:\n $(M1)\n ⋮\n $(M2)" + return nothing + end + + @test sprint(show, "text/plain", ProductManifold(M, M)) == """ + ProductManifold with 2 submanifolds: + ProductManifold($(M1), $(M2)) + ProductManifold($(M1), $(M2))""" + end + + @testset "Change Representer and Metric" begin + G = ManifoldsBase.EuclideanMetric() + @test change_representer(M, G, p1, X1) == X1 + Y1 = similar(X1) + change_representer!(M, Y1, G, p1, X1) + @test isapprox(M, p1, Y1, X1) + @test change_metric(M, G, p1, X1) == X1 + Z1 = similar(X1) + change_metric!(M, Z1, G, p1, X1) + @test isapprox(M, p1, Z1, X1) + end + + @testset "product vector transport" begin + X = log(M, p1, p2) + m = ProductVectorTransport(ParallelTransport(), ParallelTransport()) + @test default_vector_transport_method(M) === m + @test default_vector_transport_method(M, typeof(X)) === m + Y = vector_transport_to(M, p1, X, p2, m) + Y2 = similar(Y) + vector_transport_to!(M, Y2, p1, X, p2, m) + Z = -log(M, p2, p1) + @test isapprox(M, p2, Y, Z) + @test isapprox(M, p2, Y2, Z) + + Y3 = vector_transport_direction(M, p1, X, X, m) + Y4 = similar(Y) + vector_transport_direction!(M, Y4, p1, X, X, m) + Y5 = similar(Y) + vector_transport_direction!(M, Y5, p1, X, X, ParallelTransport()) + @test isapprox(M, p2, Y3, Z) + @test isapprox(M, p2, Y4, Z) + @test isapprox(M, p2, Y5, Z) + end + + @testset "Implicit product vector transport" begin + p = ArrayPartition([1, 0.0, 0.0], [4 5.0; 6 7]) + q = ArrayPartition([0.0, 1.0, 0.0], [4 8.0; 3 7.5]) + + X = log(M, p, q) + for m in [ParallelTransport(), SchildsLadderTransport(), PoleLadderTransport()] + Y = vector_transport_to(M, p, X, q, m) + Z1 = vector_transport_to( + M.manifolds[1], + submanifold_component.([p, X, q], Ref(1))..., + m, + ) + Z2 = vector_transport_to( + M.manifolds[2], + submanifold_component.([p, X, q], Ref(2))..., + m, + ) + Z = ArrayPartition(Z1, Z2) + @test isapprox(M, q, Y, Z) + Y2 = allocate(M, Y) + vector_transport_to!(M, Y2, p, X, q, m) + @test isapprox(M, q, Y2, Z) + end + for m in [ParallelTransport(), SchildsLadderTransport(), PoleLadderTransport()] + Y = vector_transport_direction(M, p, X, X, m) + Z1 = vector_transport_direction( + M.manifolds[1], + submanifold_component.([p, X, X], Ref(1))..., + m, + ) + Z2 = vector_transport_direction( + M.manifolds[2], + submanifold_component.([p, X, X], Ref(2))..., + m, + ) + Z = ArrayPartition(Z1, Z2) + @test isapprox(M, q, Y, Z) + end + end + + @testset "Parallel transport" begin + p = ArrayPartition([1, 0.0, 0.0], [4 5.0; 6 7]) + q = ArrayPartition([0.0, 1.0, 0.0], [4 8.0; 3 7.5]) + + X = log(M, p, q) + # to + Y = parallel_transport_to(M, p, X, q) + Z1 = parallel_transport_to( + M.manifolds[1], + submanifold_component.([p, X, q], Ref(1))..., + ) + Z2 = parallel_transport_to( + M.manifolds[2], + submanifold_component.([p, X, q], Ref(2))..., + ) + Z = ArrayPartition(Z1, Z2) + @test isapprox(M, q, Y, Z) + Ym = allocate(Y) + parallel_transport_to!(M, Ym, p, X, q) + @test isapprox(M, q, Y, Z) + + # direction + Y = parallel_transport_direction(M, p, X, X) + Z1 = parallel_transport_direction( + M.manifolds[1], + submanifold_component.([p, X, X], Ref(1))..., + ) + Z2 = parallel_transport_direction( + M.manifolds[2], + submanifold_component.([p, X, X], Ref(2))..., + ) + Z = ArrayPartition(Z1, Z2) + @test isapprox(M, q, Y, Z) + Ym = allocate(Y) + parallel_transport_direction!(M, Ym, p, X, X) + @test isapprox(M, q, Ym, Z) + end + + @testset "ArrayPartition" begin + @test submanifold_component(M, p1, 1) === p1.x[1] + @test submanifold_component(M, p1, Val(1)) === p1.x[1] + @test submanifold_component(p1, 1) === p1.x[1] + @test submanifold_component(p1, Val(1)) === p1.x[1] + @test submanifold_components(M, p1) === p1.x + @test submanifold_components(p1) === p1.x + end + + @testset "vee/hat" begin + X = [0.1, 0.2, 0.3, -1.0, 2.0, -3.0] + + Xc = hat(M, p1, X) + X2 = vee(M, p1, Xc) + @test isapprox(X, X2) + end + + @testset "get_coordinates" begin + # make sure `get_coordinates` does not return an `ArrayPartition` + p1 = ArrayPartition([0.0, 1.0, 0.0], [0.0 0.0; 1.0 2.0]) + X1 = ArrayPartition([1.0, 0.0, -1.0], [1.0 0.0; 2.0 1.0]) + Tp1M = TangentSpace(M, p1) + c = get_coordinates(Tp1M, p1, X1, DefaultOrthonormalBasis()) + @test c isa Vector + + p1 = ArrayPartition([0.0, 1.0, 0.0], [0.0 0.0; 3.0 20]) + X1ap = ArrayPartition([1.0, 0.0, -1.0], [1.0 0.0; 0.0 3.0]) + Tp1M = TangentSpace(M, p1) + cap = get_coordinates(Tp1M, p1, X1ap, DefaultOrthonormalBasis()) + @test cap isa Vector + end + + @testset "Basis printing" begin + p = ArrayPartition([1.0, 0.0, 0.0], [1.0 0.0; 1.0 2.0]) + B = DefaultOrthonormalBasis() + Bc = get_basis(M, p, B) + Bc_components_s = sprint.(show, "text/plain", Bc.data.parts) + @test sprint(show, "text/plain", Bc) == """ + $(typeof(B)) for a product manifold + Basis for component 1: + $(Bc_components_s[1]) + Basis for component 2: + $(Bc_components_s[2]) + """ + end + + @testset "Basis-related errors" begin + a = ArrayPartition([1.0, 0.0, 0.0], [0.0 0.0; 0.0 0.0]) + B = CachedBasis(DefaultOrthonormalBasis(), ProductBasisData(([],))) + @test_throws AssertionError get_vector!( + M, + a, + ArrayPartition([1.0, 0.0, 0.0], [0.0, 0.0]), + [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0], # this is one element too long, hence assertion error + B, + ) + @test_throws MethodError get_vector!( + M, + a, + ArrayPartition([1.0, 0.0, 0.0], [0.0 0.0; 0.0 0.0]), + [1.0, 2.0, 3.0, 4.0, 5.0, 6.0], + B, # empty elements yield a submanifold MethodError + ) + end + + @testset "Basis -- other" begin + p1 = ArrayPartition([1.0, 0.0, 0.0], [1.0 0.0; 1.0 2.0]) + X1 = ArrayPartition([0.0, 2.0, -1.0], [1.0 0.0; 2.0 1.0]) + B = DefaultOrthonormalBasis() + Bc = get_basis(M, p1, B) + + BcO = get_basis(M, p1, DiagonalizingOrthonormalBasis(X1)) + @test BcO.data isa ProductBasisData + + for basis in + [DefaultOrthonormalBasis(), get_basis(M, p1, DefaultOrthonormalBasis())] + @test length(get_vectors(M, p1, get_basis(M, p1, basis))) == 6 + X1c = get_coordinates(M, p1, X1, basis) + @test isapprox(X1c, [2.0, -1.0, 1.0, 2.0, 0.0, 1.0]) + Y1c = allocate(M, X1c) + get_coordinates!(M, Y1c, p1, X1, basis) + @test isapprox(X1c, Y1c) + @test isapprox(get_vector(M, p1, X1c, basis), X1) + Z1 = allocate(M, X1) + get_vector!(M, Z1, p1, X1c, basis) + @test isapprox(X1, X1) + end + end + + @testset "allocation promotion" begin + M2c = DefaultManifold(2; field = ℂ) + Mc = ProductManifold(M1, M2c) + @test ManifoldsBase.allocation_promotion_function(Mc, get_vector, ()) === complex + @test ManifoldsBase.allocation_promotion_function(M, get_vector, ()) === identity + end + + @testset "Riemann tensor" begin + p = ArrayPartition([0.0, 1.0, 0.0], [2.0 3.0; 0.0 0.0]) + X = ArrayPartition([1.0, 0.0, 0.0], [2.0 3.0; 0.0 0.0]) + Y = ArrayPartition([0.0, 0.0, 3.0], [-2.0 3.0; 0.0 0.0]) + Z = ArrayPartition([-1.0, 0.0, 2.0], [2.0 -3.0; 0.0 0.0]) + Xresult = ArrayPartition([6.0, 0.0, 3.0], [0.0 0.0; 0.0 0.0]) + @test isapprox(riemann_tensor(M, p, X, Y, Z), Xresult) + Xresult2 = allocate(Xresult) + riemann_tensor!(M, Xresult2, p, X, Y, Z) + @test isapprox(Xresult2, Xresult) + end + + @testset "× constructors" begin + @test M1 × M2 == M + @test M × M2 == ProductManifold(M1, M2, M2) + @test M2 × M == ProductManifold(M2, M1, M2) + @test M × M == ProductManifold(M1, M2, M1, M2) + + r1 = ExponentialRetraction() + r2 = ProjectionRetraction() + s1 = r1 × r2 + @test s1 == ProductRetraction(r1, r2) + @test "$(s1)" == "ProductRetraction($(r1), $(r2))" + @test s1 × r2 == ProductRetraction(r1, r2, r2) + @test r2 × s1 == ProductRetraction(r2, r1, r2) + @test r1 × r1 × r1 == ProductRetraction(r1, r1, r1) + @test s1 × s1 == ProductRetraction(r1, r2, r1, r2) + + ir1 = LogarithmicInverseRetraction() + ir2 = ProjectionInverseRetraction() + is1 = ir1 × ir2 + @test is1 == InverseProductRetraction(ir1, ir2) + @test "$(is1)" == "InverseProductRetraction($(ir1), $(ir2))" + @test is1 × ir2 == InverseProductRetraction(ir1, ir2, ir2) + @test ir2 × is1 == InverseProductRetraction(ir2, ir1, ir2) + @test ir1 × ir1 × ir1 == InverseProductRetraction(ir1, ir1, ir1) + @test is1 × is1 == InverseProductRetraction(ir1, ir2, ir1, ir2) + + tr1 = ParallelTransport() + tr2 = ProjectionTransport() + ts1 = tr1 × tr2 + @test ts1 == ProductVectorTransport(tr1, tr2) + @test "$(ts1)" == "ProductVectorTransport($(tr1), $(tr2))" + @test ts1 × tr2 == ProductVectorTransport(tr1, tr2, tr2) + @test tr2 × ts1 == ProductVectorTransport(tr2, tr1, tr2) + @test tr1 × tr1 × tr1 == ProductVectorTransport(tr1, tr1, tr1) + @test ts1 × ts1 == ProductVectorTransport(tr1, tr2, tr1, tr2) + end + +end diff --git a/test/runtests.jl b/test/runtests.jl index e137d91f..cf52d276 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,21 +3,19 @@ using ManifoldsBase @testset "ManifoldsBase" begin bound = 6 # six ambiguities come from possible incorrectly formed calls to `allocate` - if VERSION >= v"1.1" - num_ambiguities = length(Test.detect_ambiguities(ManifoldsBase)) - #num_ambiguities > 0 && @warn "The number of ambiguities in ManifoldsBase is $(num_ambiguities)." - if VERSION >= v"1.9-DEV" - @test num_ambiguities <= bound + 6 - elseif VERSION >= v"1.8-DEV" - @test num_ambiguities <= bound + 5 - elseif VERSION >= v"1.6-DEV" - # At the time of writing there seem to be two ambiguities regarding `getindex`, - # one with a method from SparseArrays and one from VSCode's JSON processing - # that's automatically loaded when running code in VSCode. - @test num_ambiguities <= bound + 3 - else - @test num_ambiguities == bound + 1 - end + num_ambiguities = length(Test.detect_ambiguities(ManifoldsBase)) + #num_ambiguities > 0 && @warn "The number of ambiguities in ManifoldsBase is $(num_ambiguities)." + if VERSION >= v"1.9-DEV" + @test num_ambiguities <= bound + 6 + elseif VERSION >= v"1.8-DEV" + @test num_ambiguities <= bound + 5 + elseif VERSION >= v"1.6-DEV" + # At the time of writing there seem to be two ambiguities regarding `getindex`, + # one with a method from SparseArrays and one from VSCode's JSON processing + # that's automatically loaded when running code in VSCode. + @test num_ambiguities <= bound + 3 + else + @test num_ambiguities == bound + 1 end include("decorator_traits.jl") @@ -37,4 +35,6 @@ using ManifoldsBase include("domain_errors.jl") include("vector_transport.jl") include("metric.jl") + include("product_manifold.jl") + include("fibers.jl") end diff --git a/test/test_manifolds.jl b/test/test_manifolds.jl new file mode 100644 index 00000000..571aab9e --- /dev/null +++ b/test/test_manifolds.jl @@ -0,0 +1,226 @@ +using ManifoldsBase: ℝ, ℂ, DefaultManifold, RealNumbers, EuclideanMetric +using LinearAlgebra + +# minimal sphere implementation for testing more complicated manifolds +struct TestSphere{N,𝔽} <: AbstractManifold{𝔽} end +TestSphere(N::Int, 𝔽 = ℝ) = TestSphere{N,𝔽}() + +ManifoldsBase.representation_size(::TestSphere{N}) where {N} = (N + 1,) + +function ManifoldsBase.change_representer!( + M::TestSphere, + Y, + ::ManifoldsBase.EuclideanMetric, + p, + X, +) + return copyto!(M, Y, p, X) +end + +function ManifoldsBase.change_metric!( + M::TestSphere, + Y, + ::ManifoldsBase.EuclideanMetric, + p, + X, +) + return copyto!(M, Y, p, X) +end + +function ManifoldsBase.check_point(M::TestSphere, p; kwargs...) + if !isapprox(norm(p), 1.0; kwargs...) + return DomainError( + norm(p), + "The point $(p) does not lie on the $(M) since its norm is not 1.", + ) + end + return nothing +end + +function ManifoldsBase.check_vector(M::TestSphere, p, X; kwargs...) + if !isapprox(abs(real(dot(p, X))), 0.0; kwargs...) + return DomainError( + abs(dot(p, X)), + "The vector $(X) is not a tangent vector to $(p) on $(M), since it is not orthogonal in the embedding.", + ) + end + return nothing +end + +function ManifoldsBase.exp!(M::TestSphere, q, p, X) + return exp!(M, q, p, X, one(number_eltype(X))) +end +function ManifoldsBase.exp!(::TestSphere, q, p, X, t::Number) + θ = abs(t) * norm(X) + if θ == 0 + copyto!(q, p) + else + X_scale = t * sin(θ) / θ + q .= p .* cos(θ) .+ X .* X_scale + end + return q +end + +function ManifoldsBase.get_basis_diagonalizing( + M::TestSphere{n}, + p, + B::DiagonalizingOrthonormalBasis{ℝ}, +) where {n} + A = zeros(n + 1, n + 1) + A[1, :] = transpose(p) + A[2, :] = transpose(B.frame_direction) + V = nullspace(A) + κ = ones(n) + if !iszero(B.frame_direction) + # if we have a nonzero direction for the geodesic, add it and it gets curvature zero from the tensor + V = hcat(B.frame_direction / norm(M, p, B.frame_direction), V) + κ[1] = 0 # no curvature along the geodesic direction, if x!=y + end + T = typeof(similar(B.frame_direction)) + Ξ = [convert(T, V[:, i]) for i in 1:n] + return CachedBasis(B, κ, Ξ) +end + +@inline function nzsign(z, absz = abs(z)) + psignz = z / absz + return ifelse(iszero(absz), one(psignz), psignz) +end + + +function ManifoldsBase.get_coordinates_orthonormal!(M::TestSphere, Y, p, X, ::RealNumbers) + n = manifold_dimension(M) + p1 = p[1] + cosθ = abs(p1) + λ = nzsign(p1, cosθ) + pend, Xend = view(p, 2:(n + 1)), view(X, 2:(n + 1)) + factor = λ * X[1] / (1 + cosθ) + Y .= Xend .- pend .* factor + return Y +end + +function ManifoldsBase.get_vector_orthonormal!(M::TestSphere, Y, p, X, ::RealNumbers) + n = manifold_dimension(M) + p1 = p[1] + cosθ = abs(p1) + λ = nzsign(p1, cosθ) + pend = view(p, 2:(n + 1)) + pX = dot(pend, X) + factor = pX / (1 + cosθ) + Y[1] = -λ * pX + Y[2:(n + 1)] .= X .- pend .* factor + return Y +end + +ManifoldsBase.inner(::TestSphere, p, X, Y) = dot(X, Y) + +ManifoldsBase.injectivity_radius(::TestSphere) = π + +function ManifoldsBase.inverse_retract_project!(M::TestSphere, X, p, q) + X .= q .- p + project!(M, X, p, X) + return X +end + +ManifoldsBase.is_flat(M::TestSphere) = manifold_dimension(M) == 1 + +function ManifoldsBase.log!(::TestSphere, X, p, q) + cosθ = clamp(real(dot(p, q)), -1, 1) + X .= q .- p .* cosθ + nrmX = norm(X) + if nrmX > 0 + X .*= acos(cosθ) / nrmX + end + return X +end + +ManifoldsBase.manifold_dimension(::TestSphere{N}) where {N} = N + +function ManifoldsBase.parallel_transport_to!(::TestSphere, Y, p, X, q) + m = p .+ q + mnorm2 = real(dot(m, m)) + factor = 2 * real(dot(X, q)) / mnorm2 + Y .= X .- m .* factor + return Y +end + +function ManifoldsBase.project!(::TestSphere, q, p) + q .= p ./ norm(p) + return q +end + +function ManifoldsBase.project!(::TestSphere, Y, p, X) + Y .= X .- p .* real(dot(p, X)) + return Y +end + +function Random.rand!(M::TestSphere, pX; vector_at = nothing, σ = one(eltype(pX))) + return rand!(Random.default_rng(), M, pX; vector_at = vector_at, σ = σ) +end +function Random.rand!( + rng::AbstractRNG, + M::TestSphere, + pX; + vector_at = nothing, + σ = one(eltype(pX)), +) + if vector_at === nothing + project!(M, pX, randn(rng, eltype(pX), representation_size(M))) + else + n = σ * randn(rng, eltype(pX), size(pX)) # Gaussian in embedding + project!(M, pX, vector_at, n) #project to TpM (keeps Gaussianness) + end + return pX +end + +function ManifoldsBase.riemann_tensor!(M::TestSphere, Xresult, p, X, Y, Z) + innerZX = inner(M, p, Z, X) + innerZY = inner(M, p, Z, Y) + Xresult .= innerZY .* X .- innerZX .* Y + return Xresult +end + +# from Absil, Mahony, Trumpf, 2013 https://sites.uclouvain.be/absil/2013-01/Weingarten_07PA_techrep.pdf +function ManifoldsBase.Weingarten!(::TestSphere, Y, p, X, V) + return Y .= -X * p' * V +end + +# minimal implementation of SPD (for its nontrivial metric conversion) + +struct TestSPD <: AbstractManifold{ℝ} + n::Int +end + +function ManifoldsBase.change_representer!(::TestSPD, Y, ::EuclideanMetric, p, X) + Y .= p * X * p + return Y +end + +function ManifoldsBase.change_metric!(::TestSPD, Y, ::EuclideanMetric, p, X) + Y .= p * X + return Y +end + +function ManifoldsBase.inner(::TestSPD, p, X, Y) + F = cholesky(Symmetric(convert(AbstractMatrix, p))) + return dot((F \ Symmetric(X)), (Symmetric(Y) / F)) +end + +function spd_sqrt_and_sqrt_inv(p::AbstractMatrix) + e = eigen(Symmetric(p)) + U = e.vectors + S = max.(e.values, floatmin(eltype(e.values))) + Ssqrt = Diagonal(sqrt.(S)) + SsqrtInv = Diagonal(1 ./ sqrt.(S)) + return (Symmetric(U * Ssqrt * transpose(U)), Symmetric(U * SsqrtInv * transpose(U))) +end + +function ManifoldsBase.log!(::TestSPD, X, p, q) + (p_sqrt, p_sqrt_inv) = spd_sqrt_and_sqrt_inv(p) + T = Symmetric(p_sqrt_inv * convert(AbstractMatrix, q) * p_sqrt_inv) + e2 = eigen(T) + Se = Diagonal(log.(max.(e2.values, eps()))) + pUe = p_sqrt * e2.vectors + return mul!(X, pUe, Se * transpose(pUe)) +end + +ManifoldsBase.representation_size(M::TestSPD) = (M.n, M.n) diff --git a/test/test_sphere.jl b/test/test_sphere.jl index 8572d08e..931a2e00 100644 --- a/test/test_sphere.jl +++ b/test/test_sphere.jl @@ -1,60 +1,9 @@ using LinearAlgebra using ManifoldsBase -using ManifoldsBase: ℝ, ℂ, DefaultManifold +using ManifoldsBase: ℝ, ℂ, DefaultManifold, RealNumbers using Test -# minimal sphere implementation for testing more complicated manifolds -struct TestSphere{N,𝔽} <: AbstractManifold{𝔽} end -TestSphere(N::Int, 𝔽 = ℝ) = TestSphere{N,𝔽}() - -ManifoldsBase.representation_size(::TestSphere{N}) where {N} = (N + 1,) - -function ManifoldsBase.exp!(M::TestSphere, q, p, X) - return exp!(M, q, p, X, one(number_eltype(X))) -end -function ManifoldsBase.exp!(::TestSphere, q, p, X, t::Number) - θ = norm(X) - if θ == 0 - copyto!(q, p) - else - X_scale = t * sin(θ) / θ - q .= p .* cos(θ) .+ X .* X_scale - end - return q -end - -ManifoldsBase.inner(::TestSphere, p, X, Y) = dot(X, Y) - -function ManifoldsBase.inverse_retract_project!(M::TestSphere, X, p, q) - X .= q .- p - project!(M, X, p, X) - return X -end - -function ManifoldsBase.log!(::TestSphere, X, p, q) - cosθ = clamp(real(dot(p, q)), -1, 1) - X .= q .- p .* cosθ - nrmX = norm(X) - if nrmX > 0 - X .*= acos(cosθ) / nrmX - end - return X -end - -function ManifoldsBase.project!(::TestSphere, q, p) - q .= p ./ norm(p) - return q -end - -function ManifoldsBase.project!(::TestSphere, Y, p, X) - Y .= X .- p .* real(dot(p, X)) - return Y -end - -# from Absil, Mahony, Trumpf, 2013 https://sites.uclouvain.be/absil/2013-01/Weingarten_07PA_techrep.pdf -function ManifoldsBase.Weingarten!(::TestSphere, Y, p, X, V) - return Y .= -X * p' * V -end +include("test_manifolds.jl") @testset "TestSphere" begin @testset "ShootingInverseRetraction" begin