diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index f5d739bc..dd4ba930 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -11,7 +11,7 @@ jobs: runs-on: ${{ matrix.os }} strategy: matrix: - julia-version: ["1.0", "1.6", "1.8", "~1.9.0-0"] + julia-version: ["1.0", "1.6", "1.9"] os: [ubuntu-latest, macOS-latest, windows-latest] exclude: - os: macOS-latest diff --git a/.github/workflows/documenter.yml b/.github/workflows/documenter.yml index 54abe6fe..3bcd5633 100644 --- a/.github/workflows/documenter.yml +++ b/.github/workflows/documenter.yml @@ -9,21 +9,48 @@ jobs: docs: name: Documentation runs-on: ubuntu-latest - if: "contains( github.event.pull_request.labels.*.name, 'preview docs') || github.ref == 'refs/heads/master' || contains(github.ref, 'refs/tags/')" steps: - uses: actions/checkout@v3 + - uses: quarto-dev/quarto-actions/setup@v2 + with: + version: 1.3.353 - uses: julia-actions/setup-julia@latest with: - version: 1.7 - - uses: julia-actions/julia-docdeploy@v1 + version: 1.9 + - name: Julia Cache + uses: julia-actions/cache@v1 + - name: Cache Quarto + id: cache-quarto + uses: actions/cache@v3 + env: + cache-name: cache-quarto + with: + path: tutorials/_freeze + key: ${{ runner.os }}-${{ env.cache-name }}-${{ hashFiles('tutorials/*.qmd') }} + restore-keys: | + ${{ runner.os }}-${{ env.cache-name }}- + - name: Cache Documenter + id: cache-documenter + uses: actions/cache@v3 + env: + cache-name: cache-documenter + with: + path: docs/src/tutorials + key: ${{ runner.os }}-${{ env.cache-name }}-${{ hashFiles('tutorials/*.qmd') }} + restore-keys: | + ${{ runner.os }}-${{ env.cache-name }}- + - name: Cache CondaPkg + id: cache-condaPkg + uses: actions/cache@v3 + env: + cache-name: cache-condapkg + with: + path: docs/.CondaPkg + key: ${{ runner.os }}-${{ env.cache-name }}-${{ hashFiles('docs/CondaPkg.toml') }} + restore-keys: | + ${{ runner.os }}-${{ env.cache-name }}- + - name: "Documenter rendering (including Quarto)" + run: "docs/make.jl --quarto" env: - PYTHON: "" GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} - note: - name: "Documentation deployment note." - runs-on: ubuntu-latest - if: "!contains( github.event.pull_request.labels.*.name, 'preview docs')" - steps: - - name: echo instructions - run: echo 'The Documentation is only generated and pushed on a PR if the “preview docs” label is added.' diff --git a/.github/workflows/format.yml b/.github/workflows/format.yml index 89dcab57..88abec0c 100644 --- a/.github/workflows/format.yml +++ b/.github/workflows/format.yml @@ -13,7 +13,7 @@ jobs: - uses: actions/checkout@v3 - uses: julia-actions/setup-julia@v1 with: - version: 1.4 + version: 1.9 - name: Install JuliaFormatter and format run: | julia -e 'using Pkg; Pkg.add(PackageSpec(name="JuliaFormatter"))' diff --git a/.gitignore b/.gitignore index c39388c1..c95fff7b 100644 --- a/.gitignore +++ b/.gitignore @@ -5,3 +5,6 @@ deps/deps.jl Manifest.toml docs/build +docs/src/tutorials/*.md +tutorials/_freeze +docs/.CondaPkg diff --git a/Project.toml b/Project.toml index 22a4defa..f37eb3c4 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ManifoldsBase" uuid = "3362f125-f0bb-47a3-aa74-596ffd7ef2fb" authors = ["Seth Axen ", "Mateusz Baran ", "Ronny Bergmann ", "Antoine Levitt "] -version = "0.14.6" +version = "0.14.7" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" diff --git a/bib/journal-of-the-royal-statistical-society.csl b/bib/journal-of-the-royal-statistical-society.csl new file mode 100644 index 00000000..c8d5098f --- /dev/null +++ b/bib/journal-of-the-royal-statistical-society.csl @@ -0,0 +1,305 @@ + + diff --git a/docs/CondaPkg.toml b/docs/CondaPkg.toml new file mode 100644 index 00000000..8adcf64f --- /dev/null +++ b/docs/CondaPkg.toml @@ -0,0 +1,3 @@ +[deps] +jupyter = "" +python = "3.11" \ No newline at end of file diff --git a/docs/Project.toml b/docs/Project.toml index e5a8691f..9629fb2c 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,7 +1,9 @@ [deps] +CondaPkg = "992eb4ea-22a4-4c89-a5bb-47a3300528ab" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" ManifoldsBase = "3362f125-f0bb-47a3-aa74-596ffd7ef2fb" [compat] +CondaPkg = "0.2" Documenter = "0.27" ManifoldsBase = "0.14" diff --git a/docs/make.jl b/docs/make.jl old mode 100644 new mode 100755 index fcfe26af..bdfc8d99 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,13 +1,49 @@ -using ManifoldsBase, Documenter +#!/usr/bin/env julia +# +# + +# +# (a) if docs is not the current active environment, switch to it +# (from https://github.com/JuliaIO/HDF5.jl/pull/1020/)  +if Base.active_project() != joinpath(@__DIR__, "Project.toml") + using Pkg + Pkg.activate(@__DIR__) + Pkg.develop(PackageSpec(; path = (@__DIR__) * "/../")) + Pkg.resolve() + Pkg.instantiate() +end + +# (b) Did someone say render? Then we render! +if "--quarto" ∈ ARGS + using CondaPkg + CondaPkg.withenv() do + @info "Rendering Quarto" + tutorials_folder = (@__DIR__) * "/../tutorials" + # instantiate the tutorials environment if necessary + Pkg.activate(tutorials_folder) + Pkg.resolve() + Pkg.instantiate() + Pkg.build("IJulia") # build IJulia to the right version. + Pkg.activate(@__DIR__) # but return to the docs one before + return run(`quarto render $(tutorials_folder)`) + end +end + +using Documenter: DocMeta, HTML, MathJax3, deploydocs, makedocs +using ManifoldsBase makedocs( - format = Documenter.HTML(prettyurls = false, assets = ["assets/favicon.ico"]), + format = HTML(; + mathengine = MathJax3(), + prettyurls = get(ENV, "CI", nothing) == "true", + assets = ["assets/favicon.ico"], + ), modules = [ManifoldsBase], authors = "Seth Axen, Mateusz Baran, Ronny Bergmann, and contributors.", sitename = "ManifoldsBase.jl", pages = [ "Home" => "index.md", - "How to define a manifold" => "example.md", + "How to define a manifold" => "tutorials/implement-a-manifold.md", "Design principles" => "design.md", "An abstract manifold" => "types.md", "Functions on maniolds" => [ diff --git a/docs/src/design.md b/docs/src/design.md index 45b6c772..28210e24 100644 --- a/docs/src/design.md +++ b/docs/src/design.md @@ -1,4 +1,4 @@ -# Main Design Principles +# [Main Design Principles](@id Design) The interface for a manifold is defined to be as generic as possible, such that applications can be implemented as independently as possible from an actual manifold. This way, algorithms like those from [`Manopt.jl`](https://manoptjl.org) can be implemented on _arbitrary_ manifolds. @@ -38,7 +38,7 @@ The highest layer for convenience of decorators. A usual scheme is, that a manifold might assume several things implicitly, for example the default implementation of the sphere $\mathbb S^n$ using unit vectors in $\mathbb R^{n+1}$. The embedding can be explicitly used to avoid re-implementations – the inner product can be “passed on” to its embedding. -To do so, we “decorate” the manifold by making it an [`AbstractDecoratorManifold`](@ref) and activating the right traits see [the example](@ref manifold-tutorial). +To do so, we “decorate” the manifold by making it an [`AbstractDecoratorManifold`](@ref) and activating the right traits see the tutorial [How to Implement a Manifold](@ref). The explicit case of the [`EmbeddedManifold`](@ref) can be used to distinguish different embeddings of a manifold, but also their dispatch (onto the manifold or its embedding, depending on the type of embedding) happens here. diff --git a/docs/src/example.md b/docs/src/example.md deleted file mode 100644 index 8b0c1af3..00000000 --- a/docs/src/example.md +++ /dev/null @@ -1,264 +0,0 @@ -# [How to implement your own manifold](@id manifold-tutorial) - -```@meta -CurrentModule = ManifoldsBase -``` - -## Introduction - -This tutorial explains how to implement a manifold using the `ManifoldsBase.jl` interface. -We assume that you are familiar with the basic terminology on Riemannian manifolds, especially -the dimension of a manifold, the exponential map, and the inner product on tangent spaces. -To read more about this you can for example check [^doCarmo1992], Chapter 3, first. - -Furthermore, we will look at a manifold that is isometrically embedded into a Euclidean space. - -In general you need just a data type (`struct`) that inherits from [`AbstractManifold`](@ref) to define a manifold. No function is _per se_ required to be implemented. -However, it is a good idea to provide functions that might be useful to others, for example [`check_point`](@ref check_point) and [`check_vector`](@ref check_point), as we do in this tutorial. - -In this tutorial we will - -* [model](@ref manifold-tutorial-task) the manifold -* [implement](@ref manifold-tutorial-checks) two tests, so that points and tangent vectors can be checked for validity, for example also within [`ValidationManifold`](@ref), -* [implement](@ref manifold-tutorial-fn) two functions, the exponential map and the manifold dimension. -* [decorate](@ref manifold-tutorial-emb) the manifold with an embedding to gain further features. - -The next two sections are just a technical detail and the necessary `import`s to extend the functions defined in this interface. - -## [Technical preliminaries](@id manifold-tutorial-prel) - -There are only two small technical things we need to explain at this point before we get started. -First of all our [`AbstractManifold`](@ref)`{𝔽}` has a parameter `𝔽`. -This parameter indicates the [`number_system`](@ref) the manifold is based on, for example `ℝ` for real manifolds, which is short for [`RealNumbers`](@ref)`()` or `ℂ` for complex manifolds, a shorthand for [`ComplexNumbers`](@ref)`()`. - -Second, this interface usually provides both an allocating and an in-place variant of each function, for example for the [`exp`](@ref)onential map [implemented below](@ref manifold-tutorial-fn) this interface provides `exp(M, p, X)` to compute the exponential map and `exp!(M, q, p, X)` to compute the exponential map in the memory provided by `q`, mutating that input. -the convention is, that the manifold is the first argument -- in both function variants -- the in-place variant then has the input to be mutated in second place, and the remaining parameters are again the same (`p`and `X` here). -We usually refer to these two variants of the same function as the allocating (`exp`) function and the in-place (`exp!`) one. - -The convention for this interface is to __document the allocation function__, which by default allocates the necessary memory and calls the in-place function. So the convention is to just __implement the in-place function__, unless there is a good reason to provide an implementation for both. -For more details see [the design section on in-place and non-mutating functions](@ref inplace-and-noninplace) - -For performance reasons scaled variants of retractions `retraction!(M, q, p, X, t, m)` and exponential maps `exp!(M, q, p, X, t)` should be implemented. Scale is specified by an optional argument that comes after the tangent vector and by default it is equal to 1. The variant without scaling, e.g. `exp!(M, q, p, X)`, can be implemented as well if there is a good reason like performance, but the default fallback of this variant is to call the previous one with `t=1`. - -## [Startup](@id manifold-tutorial-startup) - -As a start, let's load `ManifoldsBase.jl` and import the functions we consider throughout this tutorial. - -```@example manifold-tutorial -using ManifoldsBase, LinearAlgebra, Test -import ManifoldsBase: check_point, check_vector, manifold_dimension, exp!, inner, representation_size, get_embedding -import Base: show -``` - -We load `LinearAlgebra` for some computations. `Test` is only loaded for illustrations in the examples. - -We import the in-place variant of the [`exp`](@ref)onential map, as just discussed above. - -## [The manifold](@id manifold-tutorial-task) - -The manifold we want to implement here a sphere, with radius $r$. -Since this radius is a property inherent to the manifold, it will become a field of the manifolds `struct`. -The second information, we want to store is the dimension of the sphere, for example whether it's the 1-sphere, i.e. the circle, represented by vectors $p\in\mathbb R^2$ of norm $r$ or the 2-sphere in $\mathbb R^3$ of radius $r$. -Since the latter might be something we want to [dispatch](https://en.wikipedia.org/wiki/Multiple_dispatch) on, we model it as a parameter of the type. -In general the `struct` of a manifold should provide information about the manifold, which are inherent to the manifold or has to be available without a specific point or tangent vector present. -This is -- most prominently -- all information required to determine the manifold dimension. - -Note that this a slightly more general manifold than the [Sphere](https://juliamanifolds.github.io/Manifolds.jl/stable/manifolds/sphere.html) in [Manifolds.jl](https://juliamanifolds.github.io/Manifolds.jl/stable/index.html) - -For our example we define the following `struct`. -While a first implementation might also just take [`AbstractManifold`](@ref)`{ℝ}` as supertype, we directly take -[`AbstractDecoratorManifold`](@ref)`{ℝ}`, which will be useful later on. -For now it does not make a difference. - -```@example manifold-tutorial -""" - ScaledSphere{N} <: AbstractDecoratorManifold{ℝ} - -Define an `N`-sphere of radius `r`. Construct by `ScaledSphere(radius,n)`. -""" -struct ScaledSphere{N} <: AbstractDecoratorManifold{ManifoldsBase.ℝ} where {N} - radius::Float64 -end -ScaledSphere(radius, n) = ScaledSphere{n}(radius) -Base.show(io::IO, M::ScaledSphere{n}) where {n} = print(io, "ScaledSphere($(M.radius),$n)") -nothing #hide -``` - -Here, the last line just provides a nicer print of a variable of that type. -Now we can already initialize our manifold that we will use later, the $2$-sphere of radius $1.5$. - -```@example manifold-tutorial -S = ScaledSphere(1.5, 2) -``` - -## [Checking points and tangents](@id manifold-tutorial-checks) - -Points on a manifold are usually represented as vector, matrices or more generally arrays. -Since we consider vectors of a certain norm (and space dimension), our points are vectors. -For an arbitrary vector we would first like to check, that it is a valid point on the manifold. -For this one can use the function [`is_point`](@ref is_point(M::AbstractManifold, p; kwargs...)). -This is a function on [layer 1](@ref design-layer1) which handles special cases as well cases, so it should not be implemented directly by a user of this interface. -The functions that have to be implemented can be found on [layer 3](@ref design-layer3). Generically, for both [`is_point`](@ref is_point(M::AbstractManifold, p; kwargs...)) and [`is_vector`](@ref is_vector(M::AbstractManifold, p, X; kwargs...)), this layer contains a function to check correct size of an array, called [`check_size`](@ref ManifoldsBase.check_size) -For the test of points the function to implement is [`check_point`](@ref ManifoldsBase.check_point) which we actually will implement, analogously there exists also [`check_vector`](@ref ManifoldsBase.check_vector). -These functions return `nothing` if the point (vector, size) is a correct/valid and returns an error (but not throw it) otherwise. -This is usually a `DomainError`. - -We have to check two things: that a point `p` is a vector with `N+1` entries and its norm is the desired radius. - -A first thing we have specify is how points and tangent vectors are represented, that is we have to specify their [`representation_size`](@ref) - -```@example manifold-tutorial -representation_size(::ScaledSphere{N}) where {N} = N+1 -nothing #hide -``` - -This already finishes the size check which [`check_size`](@ref ManifoldsBase.check_size) performs by default (based on the representation size). - -If something has to only hold up to precision, we can pass that down, too using the `kwargs...`, so all three `check_` functions should usually have these in their signature. -For our manifold we have to check that the norm of a point `p` is approximately the specified `radius`. - -```@example manifold-tutorial -function check_point(M::ScaledSphere{N}, p; kwargs...) where {N} - if !isapprox(norm(p), M.radius; kwargs...) - return DomainError(norm(p), "The norm of $p is not $(M.radius).") - end - return nothing -end -nothing #hide -``` - -Similarly, we can verify, whether a tangent vector `X` is valid. -Its size is again already checked using [`check_size`](@ref ManifoldsBase.check_size), -so the only remaining property to verify is, that `X` is orthogonal to `p`. -We can again use the `kwargs`. - -```@example manifold-tutorial -function check_vector(M::ScaledSphere, p, X; kwargs...) - if !isapprox(dot(p,X), 0.0; kwargs...) - return DomainError(dot(p,X), "The tangent $X is not orthogonal to $p.") - end - return nothing -end -nothing #hide -``` - -Note that the function [`is_vector`](@ref is_vector(M::AbstractManifold, p, X; kwargs...)) -even can check that the base point of `X` (the `p` the tangent space belongs to), can be checked for validity, -see its keyword argument `check_base_point`, so within [`check_vector`](@ref ManifoldsBase.check_vector) -this can be (implicitly) assumed to hold. - -to test points we can now use - -```@example manifold-tutorial -is_point(S, [1.0,0.0,0.0]) # norm 1, so not on S, returns false -@test_throws DomainError is_point(S, [1.5,0.0], true) # only on R^2, throws an error. -p = [1.5,0.0,0.0] -X = [0.0,1.0,0.0] -# The following two tests return true -[ is_point(S, p); is_vector(S,p,X) ] -``` - -## [Functions on the manifold](@id manifold-tutorial-fn) - -For the [`manifold_dimension`](@ref manifold_dimension(M::AbstractManifold)) we have to just return the `N` parameter - -```@example manifold-tutorial -manifold_dimension(::ScaledSphere{N}) where {N} = N -manifold_dimension(S) -``` - -Note that we can even omit the variable name in the first line since we do not have to access any field or use the variable otherwise. - -To implement the [`exp`](@ref)onential map, we have to implement the formula for great arcs, given a start point `p` and a direction `X` on the $n$-sphere of radius $r$ the formula reads - -````math -\exp_p X = \cos\Bigl(\frac{1}{r}\lVert X \rVert\Bigr)p + \sin\Bigl(\frac{1}{r}\lVert X \rVert\Bigr)\frac{r}{\lVert X \rVert}X. -```` - -Note that with this choice we for example implicitly assume that the manifold is equipped with that certain metric. -This is the default within this interface. -To distinguish different metrics, see [`MetricManifold`](https://juliamanifolds.github.io/Manifolds.jl/latest/manifolds/metric.html) in [Manifolds.jl](https://juliamanifolds.github.io/Manifolds.jl/stable/) for more details. -Since we here only consider one metric, we do not have to specify that. - -An implementation of the mutation version, see the [technical note](@ref manifold-tutorial-prel) for the naming and reasoning, reads - -```@example manifold-tutorial -function exp!(M::ScaledSphere{N}, q, p, X, t::Number) where {N} - nX = abs(t) * norm(X) - if nX == 0 - q .= p - else - q .= cos(nX/M.radius)*p + M.radius*sin(nX/M.radius) .* t .* (X./nX) - end - return q -end -function exp!(M::ScaledSphere{N}, q, p, X) where {N} - exp!(M, q, p, X, 1) -end -nothing #hide -``` - -Two variants are implemented above: one with the scaling argument `t` and one without. It isn't always necessary to implement both but doing so reduces the chance of ambiguity errors. - -A first easy check can be done taking `p` from above and any vector `X` of length `1.5π` from its tangent space. -The resulting point is opposite of `p`, i.e. `-p` and it is of course a valid point on `S`. - -```@example manifold-tutorial -q = exp(S, p, [0.0,1.5π,0.0]) -[isapprox(p, -q); is_point(S, q)] -``` - -## [Adding an isometric embedding](@id manifold-tutorial-emb) - -Since the sphere is isometrically embedded, we do not have to implement the [`inner`](@ref)`(M,p,X,Y)` for tangent vectors `X`, `Y` in the tangent space at `p` , but we can “delegate” it to the embedding. The embedding is the [Euclidean](https://juliamanifolds.github.io/Manifolds.jl/latest/manifolds/euclidean.html). -The same manifold with a little smaller feature set is available in `ManifoldsBase.jl` as `DefaultManifold` for testing purposes. - -```@example manifold-tutorial -using ManifoldsBase: DefaultManifold, IsIsometricEmbeddedManifold -import ManifoldsBase: active_traits, get_embedding -using ManifoldsBase: merge_traits -``` - -Now we can activate a decorator by specifying that the sphere has the [`IsIsometricEmbeddedManifold`](@ref) trait for the functions `f` on our scaled sphere manifold by writing - -```@example manifold-tutorial -active_traits(f, ::ScaledSphere, args...) = merge_traits(IsIsometricEmbeddedManifold()) -nothing #hide -``` - -and then specifying that said embedding is the `DefaultManifold`. - -```@example manifold-tutorial -get_embedding(::ScaledSphere{N}) where {N} = DefaultManifold(N+1) -nothing #hide -``` - -Now metric related functions are passed to this embedding, for example the inner product. -It now works by using the inner product from the embedding, so we can compute the inner product by calling [`inner`](@ref) - -```@example manifold-tutorial -X = [0.0, 0.1, 3.0] -Y = [0.0, 4.0, 0.2] -# returns 1.0 by calling the inner product in DefaultManifold(3) -inner(S, p, X, Y) -``` - -## [Conclusion](@id manifold-tutorial-outlook) - -You can now just continue implementing further functions from `ManifoldsBase.jl` -but with just [`exp!`](@ref exp!(M::AbstractManifold, q, p, X, t::Number)) you for example already have - -* [`geodesic`](@ref geodesic(M::AbstractManifold, p, X)) the (not necessarily shortest) geodesic emanating from `p` in direction `X`. -* the [`ExponentialRetraction`](@ref), that the [`retract`](@ref retract(M::AbstractManifold, p, X)) function uses by default. - -For the [`shortest_geodesic`](@ref shortest_geodesic(M::AbstractManifold, p, q)) the implementation of a logarithm [`log`](@ref ManifoldsBase.log(M::AbstractManifold, p, q)), or just [`log!`](@ref log!(M::AbstractManifold, X, p, q)) is sufficient. - -Sometimes a default implementation is provided; for example if you implemented [`inner`](@ref inner(M::AbstractManifold, p, X, Y)), the [`norm`](@ref norm(M, p, X)) is defined. You should overwrite it, if you can provide a more efficient version. For a start the default should suffice. -With [`log!`](@ref log!(M::AbstractManifold, X, p, q)) and [`inner`](@ref inner(M::AbstractManifold, p, X, Y)) you get the [`distance`](@ref distance(M::AbstractManifold, p, q)), and so. - -In summary with just these few functions you can already explore the first things on your own manifold. Whenever a function from `Manifolds.jl` requires another function to be specifically implemented, you get a reasonable error message. - -## Literature - -[^doCarmo1992]: - > do Carmo, Manfredo __Riemannian Geometry__, Birkhäuser Boston, 1992, ISBN: 0-8176-3490-8. diff --git a/src/metric.jl b/src/metric.jl index cbdace04..508c9386 100644 --- a/src/metric.jl +++ b/src/metric.jl @@ -29,14 +29,14 @@ abstract type RiemannianMetric <: AbstractMetric end EuclideanMetric <: RiemannianMetric A general type for any manifold that employs the Euclidean Metric, for example -the [`Euclidean`](@ref) manifold itself, or the [`Sphere`](@ref), where every +the [`Euclidean`](https://juliamanifolds.github.io/Manifolds.jl/latest/manifolds/euclidean.html) manifold itself, or the [`Sphere`](https://juliamanifolds.github.io/Manifolds.jl/latest/manifolds/sphere.html), where every tangent space (as a plane in the embedding) uses this metric (in the embedding). Since the metric is independent of the field type, this metric is also used for the Hermitian metrics, i.e. metrics that are analogous to the `EuclideanMetric` but where the field type of the manifold is `ℂ`. -This metric is the default metric for example for the [`Euclidean`](@ref) manifold. +This metric is the default metric for example for the [`Euclidean`](https://juliamanifolds.github.io/Manifolds.jl/latest/manifolds/euclidean.html) manifold. """ struct EuclideanMetric <: RiemannianMetric end diff --git a/src/retractions.jl b/src/retractions.jl index a6a2cdb5..4892db13 100644 --- a/src/retractions.jl +++ b/src/retractions.jl @@ -21,7 +21,7 @@ An abstract type for representing approximate inverse retraction methods. abstract type ApproximateInverseRetraction <: AbstractInverseRetractionMethod end """ - ApproximateRetraction <: AbstractInverseRetractionMethod + ApproximateRetraction <: AbstractRetractionMethod An abstract type for representing approximate retraction methods. """ @@ -110,7 +110,12 @@ end PolarRetraction <: AbstractRetractionMethod Retractions that are based on singular value decompositions of the matrix / matrices -for point and tangent vector on a [`AbstractManifold`](@ref) +for point and tangent vectors. + +!!! note "Technical Note" + Though you would call e.g. [`retract`](@ref)`(M, p, X, PolarRetraction())`, + to implement a polar retraction, define [`retract_polar`](@ref)`(M, p, X, t)` + or [`retract_polar!`](@ref)`(M, q, p, X, t)` for your manifold `M`. """ struct PolarRetraction <: AbstractRetractionMethod end @@ -118,6 +123,11 @@ struct PolarRetraction <: AbstractRetractionMethod end ProjectionRetraction <: AbstractRetractionMethod Retractions that are based on projection and usually addition in the embedding. + +!!! note "Technical Note" + Though you would call e.g. [`retract`](@ref)`(M, p, X, ProjectionRetraction())`, + to implement a projection retraction, define [`retract_project`](@ref)`(M, p, X, t)` + or [`retract_project!`](@ref)`(M, q, p, X, t)` for your manifold `M`. """ struct ProjectionRetraction <: AbstractRetractionMethod end @@ -126,6 +136,11 @@ struct ProjectionRetraction <: AbstractRetractionMethod end Retractions that are based on a QR decomposition of the matrix / matrices for point and tangent vector on a [`AbstractManifold`](@ref) + +!!! note "Technical Note" + Though you would call e.g. [`retract`](@ref)`(M, p, X, QRRetraction())`, + to implement a QR retraction, define [`retract_qr`](@ref)`(M, p, X, t)` + or [`retract_qr!`](@ref)`(M, q, p, X, t)` for your manifold `M`. """ struct QRRetraction <: AbstractRetractionMethod end @@ -162,6 +177,11 @@ end SoftmaxRetraction <: AbstractRetractionMethod Describes a retraction that is based on the softmax function. + +!!! note "Technical Note" + Though you would call e.g. [`retract`](@ref)`(M, p, X, SoftmaxRetraction())`, + to implement a softmax retraction, define [`retract_softmax`](@ref)`(M, p, X, t)` + or [`retract_softmax!`](@ref)`(M, q, p, X, t)` for your manifold `M`. """ struct SoftmaxRetraction <: AbstractRetractionMethod end @@ -169,7 +189,15 @@ struct SoftmaxRetraction <: AbstractRetractionMethod end @doc raw""" PadeRetraction{m} <: AbstractRetractionMethod -A retraction based on the Padé approximation of order $m$ +A retraction based on the Padé approximation of order ``m`` + +# Constructor + PadeRetraction(m::Int) + +!!! note "Technical Note" + Though you would call e.g. [`retract`](@ref)`(M, p, X, PadeRetraction(m))`, + to implement a Padé retraction, define [`retract_pade`](@ref)`(M, p, X, t, m)` + or [`retract_pade!`](@ref)`(M, q, p, X, t, m)` for your manifold `M`. """ struct PadeRetraction{m} <: AbstractRetractionMethod end @@ -184,6 +212,12 @@ end A retraction based on the Cayley transform, which is realized by using the [`PadeRetraction`](@ref)`{1}`. + +!!! note "Technical Note" + Though you would call e.g. [`retract`](@ref)`(M, p, X, CayleyRetraction())`, + to implement a caley retraction, define [`retract_cayley`](@ref)`(M, p, X, t)` + or [`retract_cayley!`](@ref)`(M, q, p, X, t)` for your manifold `M`. + By default both these functions fall back to calling a [`PadeRetraction`](@ref)`(1)`. """ const CayleyRetraction = PadeRetraction{1} @@ -213,9 +247,14 @@ struct LogarithmicInverseRetraction <: AbstractInverseRetractionMethod end @doc raw""" - PadeInverseRetraction{m} <: AbstractRetractionMethod + PadeInverseRetraction{m} <: AbstractInverseRetractionMethod 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))`, + to implement an inverse Padé retraction, define [`inverse_retract_pade`](@ref)`(M, p, q, m)` + or [`inverse_retract_pade!`](@ref)`(M, X, p, q, m)` for your manifold `M`. """ struct PadeInverseRetraction{m} <: AbstractInverseRetractionMethod end @@ -230,6 +269,12 @@ end A retraction based on the Cayley transform, which is realized by using the [`PadeRetraction`](@ref)`{1}`. + +!!! note "Technical Note" + Though you would call e.g. [`inverse_retract`](@ref)`(M, p, q, CayleyInverseRetraction())`, + to implement an inverse caley retraction, define [`inverse_retract_cayley`](@ref)`(M, p, q)` + or [`inverse_retract_cayley!`](@ref)`(M, X, p, q)` for your manifold `M`. + By default both these functions fall back to calling a [`PadeInverseRetraction`](@ref)`(1)`. """ const CayleyInverseRetraction = PadeInverseRetraction{1} @@ -238,6 +283,11 @@ const CayleyInverseRetraction = PadeInverseRetraction{1} Inverse retractions that are based on a singular value decomposition of the matrix / matrices for point and tangent vector on a [`AbstractManifold`](@ref) + +!!! note "Technical Note" + Though you would call e.g. [`inverse_retract`](@ref)`(M, p, q, PolarInverseRetraction())`, + to implement an inverse polar retraction, define [`inverse_retract_polar`](@ref)`(M, p, q)` + or [`inverse_retract_polar!`](@ref)`(M, X, p, q)` for your manifold `M`. """ struct PolarInverseRetraction <: AbstractInverseRetractionMethod end @@ -245,6 +295,11 @@ struct PolarInverseRetraction <: AbstractInverseRetractionMethod end ProjectionInverseRetraction <: AbstractInverseRetractionMethod Inverse retractions that are based on a projection (or its inversion). + +!!! note "Technical Note" + Though you would call e.g. [`inverse_retract`](@ref)`(M, p, q, ProjectionInverseRetraction())`, + to implement an inverse projection retraction, define [`inverse_retract_project`](@ref)`(M, p, q)` + or [`inverse_retract_project!`](@ref)`(M, X, p, q)` for your manifold `M`. """ struct ProjectionInverseRetraction <: AbstractInverseRetractionMethod end @@ -253,6 +308,11 @@ struct ProjectionInverseRetraction <: AbstractInverseRetractionMethod end Inverse retractions that are based on a QR decomposition of the matrix / matrices for point and tangent vector on a [`AbstractManifold`](@ref) + +!!! note "Technical Note" + Though you would call e.g. [`inverse_retract`](@ref)`(M, p, q, QRInverseRetraction())`, + to implement an inverse QR retraction, define [`inverse_retract_qr`](@ref)`(M, p, q)` + or [`inverse_retract_qr!`](@ref)`(M, X, p, q)` for your manifold `M`. """ struct QRInverseRetraction <: AbstractInverseRetractionMethod end @@ -306,7 +366,7 @@ end """ - InverseRetractionWithKeywords{R<:AbstractRetractionMethod,K} <: AbstractRetractionMethod + InverseRetractionWithKeywords{R<:AbstractRetractionMethod,K} <: AbstractInverseRetractionMethod Since inverse retractions might have keywords, this type is a way to set them as an own type to be used as a specific inverse retraction. @@ -342,6 +402,11 @@ end SoftmaxInverseRetraction <: AbstractInverseRetractionMethod Describes an inverse retraction that is based on the softmax function. + +!!! note "Technical Note" + Though you would call e.g. [`inverse_retract`](@ref)`(M, p, q, SoftmaxInverseRetraction())`, + to implement an inverse softmax retraction, define [`inverse_retract_softmax`](@ref)`(M, p, q)` + or [`inverse_retract_softmax!`](@ref)`(M, X, p, q)` for your manifold `M`. """ struct SoftmaxInverseRetraction <: AbstractInverseRetractionMethod end diff --git a/tutorials/.gitignore b/tutorials/.gitignore new file mode 100644 index 00000000..075b2542 --- /dev/null +++ b/tutorials/.gitignore @@ -0,0 +1 @@ +/.quarto/ diff --git a/tutorials/Project.toml b/tutorials/Project.toml new file mode 100644 index 00000000..53127daa --- /dev/null +++ b/tutorials/Project.toml @@ -0,0 +1,11 @@ +[deps] +IJulia = "7073ff75-c697-5162-941a-fcdaad2a7d2a" +Manifolds = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e" +ManifoldsBase = "3362f125-f0bb-47a3-aa74-596ffd7ef2fb" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" + +[compat] +IJulia = "1" +Manifolds = "0.8.46" +ManifoldsBase = "0.14.5" +Plots = "1.38" diff --git a/tutorials/_quarto.yml b/tutorials/_quarto.yml new file mode 100644 index 00000000..50bbc704 --- /dev/null +++ b/tutorials/_quarto.yml @@ -0,0 +1,27 @@ +project: + title: "ManifoldsBase.jl Tutorials" + output-dir: ../docs/src/tutorials + render: + - "*.qmd" + +crossref: + fig-prefix: Figure + tbl-prefix: Table +bibliography: + - ../CITATION.bib +csl: ../bib/journal-of-the-royal-statistical-society.csl +fig-format: png + +execute: + freeze: auto + eval: true + echo: true + output: true + + +format: + commonmark: + variant: -raw_html+tex_math_dollars + wrap: preserve + +jupyter: julia-1.9 \ No newline at end of file diff --git a/tutorials/implement-a-manifold.qmd b/tutorials/implement-a-manifold.qmd new file mode 100644 index 00000000..5da80cd1 --- /dev/null +++ b/tutorials/implement-a-manifold.qmd @@ -0,0 +1,283 @@ +--- +title: "How to Implement a Manifold" +--- + +````{=commonmark} +```@meta +CurrentModule = ManifoldsBase +``` +```` + +This tutorial illustrates, how to implement your very first manifold. +We start from the very beginning and cover the basic ideas of the interface provided by +`ManifoldsBase.jl` interface. + +```{julia} +#| echo: false +#| code-fold: true +#| output: false +using Pkg; +cd(@__DIR__) +Pkg.activate("."); # for reproducibility use the local tutorial environment. +``` + +## Preliminaries + +We will use a simple example in this tutorial, since the main focus here is to illustrate how to define a manifold. +We will use the sphere of radius $r$ embedded in $\mathbb R^{d+1}$, i.e. all vectors of length $r$. Formally we define + +```math +\mathbb S_r^d := +\bigl\{ + p \in \mathbb R^{d+1} + \big| + \lVert p \rVert = r +\bigr\} +``` + +When defining a Riemannian manifold mathematically, there is several things to keep in mind, for example the metric imposed on the tangent spaces. For this interface we assume these things to be given implicitly for a first implementation, but they can be made more precise when necessary. + +The only thing we have to be aware of for now is the [`number_system`](@ref), i.e. whether our manifold is a real-valued or a complex-valued manifold. +The abstract type all manifolds inherit from, the [`AbstractManifold`](@ref)`{𝔽}` has this number system as a parameter. +The usual parameter we will use are the [`RealNumbers`](@ref)`()`, +which have a short hand in `ManifoldsBase.jl`, namely `ℝ`. +The second one are the [`ComplexNumbers`](@ref)`()`, or `ℂ` for short. + +```{julia} +using LinearAlgebra, ManifoldsBase +using ManifoldsBase: ℝ +``` + +## Defining a manifold + +A manifold itself is a `struct` that is a subtype of `AbstractManifold` and should contain. We usually recommend to also document your new manifold. + +Since the [`Sphere`](https://juliamanifolds.github.io/Manifolds.jl/latest/manifolds/sphere.html) is already a name used within [`Manifolds.jl`](https://github.com/JuliaManifolds/Manifolds.jl/), let's use a slightly more specific name. We define + +```{julia} +#| output: false +""" + ScaledSphere <: AbstractManifold{ℝ} + +Define a sphere of fixed radius + +# Fields + +* `dimension` dimension of the sphere +* `radius` the radius of the sphere + +# Constructor + + ScaledSphere(dimension,radius=1.0) + +Initialize the manifold to a certain `dimension` and `radius`, +which by default is set to `1.0` +""" +struct ScaledSphere <: AbstractManifold{ℝ} + dimension::Int + radius::Float64 +end +``` + +And we can directly use this manifold and set + +```{julia} +M = ScaledSphere(2,1.5) +``` + +## Functions I: Manifold properties + +While the interface provides a lot of possible functions to define for your manifold, you only need to define those that are necessary for your implementation. +If you are using other packages depending on `ManifoldsBase.jl` you will often just get a “Method not defined” and sometimes an ambiguity error indicating that a function is missing that is required for a certain task. + +We can first start with a technical function which internally ist often used. Any of our points or tangent vectors is represented as a $(d+1)$-dimensional vector. This is internally often used when allocating memory, see [`representation_size`](@ref). +It returns a tuple representing the size of arrays for valid points. +We define + +```{julia} +#| output: false +import ManifoldsBase: representation_size +representation_size(M::ScaledSphere) = (M.dimension+1,) +``` + +Similarly, we can implement the function returning the dimension of the manifold, cf. [`manifold_dimension`](@ref) as + +```{julia} +#| output: false +import ManifoldsBase: manifold_dimension +manifold_dimension(M::ScaledSphere) = M.dimension +``` + +and we can now easily use them to access the dimension of the manifold + +```{julia} +manifold_dimension(M) +``` + +## Functions II: Verifying Points and tangent vectors + + +The first two functions we want to define are those to check points and tangent vectors for our manifold. +Let's first clarify what the tangent space looks like. The directions “we can walk into” from a point $p\in \mathbb S_r^d$ are all $X$ that are orthogonal to $p$, which is the plane/vector space tangent to the sphere. Formally + +```math +T_p\mathbb S_r^d := +\bigl\{ + X \in \mathbb R^{d+1} + \big| + \langle p, X \rangle = 0 +\bigr\}, \qquad p \in \mathbb S_r^d +``` + +to verify either `p` or `X` one uses [`is_point`](@ref)`(M,p)` +and [`is_vector`](@ref)`(M, p, X)` respectively. Since both involve some automatic options and possibililities, for example whether to throw an error or just return false, both mention that the actual functions to implement are [`check_point`](@ref) and [`check_vector`](@ref), which both do not throw but _return_ an error if something is wrong. + +In principle we would have to check two properties, namely that the size of `p` is of correct size `M.dimension+1` and that its norm is `M.radius`. Luckily, by defining [`representation_size`](@ref) the first check is automatically done already – actually for both points and vectors. +We define + +```{julia} +#| output: false +import ManifoldsBase: check_point +function check_point(M::ScaledSphere, p; kwargs...) + if !isapprox(norm(p), M.radius; kwargs...) + return DomainError(norm(p), "The norm of $p is not $(M.radius).") + end + return nothing +end +``` + +And we can directly test the function. To see all 3 failing ones, we switch from errors to warnings in the check + +```{julia} +is_point(M, [1.5, 0.0], :warn) # wrong size +is_point(M, [1.0, 0.0, 0.0], :warn) # wrong norm +is_point(M, [1.5, 0.0, 0.0], :warn) # on the manifold, returns true +``` + +similarly for vectors, we just have to implement the orthogonality check. + +```{julia} +#| output: false +import ManifoldsBase: check_vector +function check_vector(M::ScaledSphere, p, X; kwargs...) + if !isapprox(dot(p,X), 0.0; kwargs...) + return DomainError( + dot(p,X), + "The tangent vector $X is not orthogonal to $p." + ) + end + return nothing +end +``` + +and again, the high level interface can be used to display warning for vectors not fulfilling the criterion, but we can also +activate a check for the point using the last positional argument + +```{julia} +is_vector(M, [1.5, 0.0, 0.0], [0.0, 1.0], :warn) # wrong size +is_vector(M, [1.5, 0.0, 0.0], [1.0, 1.0, 0.0], :warn) # not orthogonal norm +is_vector(M, [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], :warn, true) # point not valid +is_vector(M, [1.5, 0.0, 0.0], [0.0, 1.0, 0.0], :warn, true) # returns true +``` + + +## Functions on Manifolds III: The exponential map and a retraction. + +For the final group of functions, we want to implement the [`exp`](@ref)onential map and a [`retract`](@ref)ion. Both are ways to “move around” on the manifold, that is, given a point $p$ and a +tangent vector indicating a “walking direction”, the two functions we define will return a new point $q$ on the manifold. + +For functions that compute a new point or tangent vector, `ManifoldsBase.jl` always provides two varinats: One that allocates new memory and one, that allows to provide the memory, the result should be returned in – to spare memory allocations. + +Let's first take a look at what the exponential map is defined like. We follow the shortest curves, that is great arcs, on the sphere. Formally we have + +````math +\exp_p X = +\cos\Bigl(\frac{1}{r}\lVert X \rVert\Bigr)p + +\sin\Bigl(\frac{1}{r}\lVert X \rVert\Bigr)\frac{X}{\lVert X \rVert}. +```` + +In fact, from the two functions above, [`exp`](@ref)`(M, p, X)` and [`exp!`](@ref)`(M, q, p, X)`, that works in place of `q`, we only have to implement the second. +The first one, `exp` by default falls back to allocating memory and calling the secnod. Sp `exp` should only be defined, if there is a special reason for. +Furthermore, we usually do not verify/check inputs to spare time. If a user feels insecure, they could for example use the [`ValidationManifold`](@ref) wrapper which adds explicitly checks of inputs and outputs. + +We define + +```{julia} +#| output: false +import ManifoldsBase: exp! +function exp!(M::ScaledSphere, q, p, X) + nX = norm(X) + if nX == 0 + q .= p + else + q .= cos(nX/M.radius)*p + M.radius*sin(nX/M.radius) .* (1/nX) .* X + end + return q +end +``` + +and we can directly test our function starting in the north pole and “waling down” to the equator + +```{julia} +q = exp(M, [0.0, 0.0, 1.5], [0.75π, 0.0, 0.0]) +``` + +but we also get the other variants for free, for example + +```{julia} +q2 = zero(q) +exp!(M, q2, [0.0, 0.0, 1.5], [0.75π, 0.0, 0.0]) +q2 +``` + +or the one that shortens or enlongates the path by a factor, +for example, if we walk twice the distance, we reach the opposite point + +```{julia} +exp!(M, q2, [0.0, 0.0, 1.5], [0.75π, 0.0, 0.0], 2.0) +q2 +``` + +Of course we can easliy check that both points we computed still lie on the sphere + +```{julia} +all([is_point(M, q), is_point(M, q2)]) +``` + +Since the exponential map might in general be expensive, we can do a similar implementation with the [`ProjectionRetraction`](@ref). +Here, we really have to take into account, that the interface is ``[designed with 3 levels](@ref Design)``{=commonmark} in mind: +While the actual function we would call in the end is `retract(M, p, X, ProjectionRetraction())` (or its `!` variant), we actually have to implement `retract_project!(M, q, p, X, t)` for technical details, that are a bit beyond this introductionary tutorial. In short this split avoids ambiguity errors for decorators of the manifolds. We define + +```{julia} +import ManifoldsBase: retract_project! +function retract_project!(M::ScaledSphere, q, p, X, t) + q .= (p + t*X) .* (M.radius/norm(p + t*X)) + return q +end +``` + +And to test also this function, and avoiding allocations at the same time, we call + +```{julia} +retract!(M, q, [0.0, 0.0, 1.5], [0.75π, 0.0, 0.0], ProjectionRetraction()) +``` + +Finally, there is [`default_retraction_method`](@ref) to specify which is the default retraction to use. By default this is + +```{julia} +default_retraction_method(M) +``` + +But we can easily specify this for our manifold as well, for example defining + +```{julia} +import ManifoldsBase: default_retraction_method +default_retraction_method(::ScaledSphere) = ProjectionRetraction() +``` + +Then +```{julia} +default_retraction_method(M) +``` +and retract without a method specified would always fall back to using the projection retraction instead of the exponential map. +Note that for compatibilty there is the [`AbstractRetractionMethod`](@ref) called [`ExponentialRetraction`](@ref) which makes [`retract`](@ref) fall back to calling [`exp`](@ref).