From 58f6a6b00a5c44bc940d742efbca01a4b9856e9b Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Sat, 16 Mar 2024 09:47:41 +0100 Subject: [PATCH 01/14] Sectional curvature --- NEWS.md | 5 ++ Project.toml | 4 +- docs/make.jl | 1 + docs/src/features/utilities.md | 12 +-- docs/src/misc/notation.md | 1 + docs/src/references.bib | 10 +++ src/Manifolds.jl | 7 ++ src/manifolds/Euclidean.jl | 27 +++++++ src/manifolds/Hyperbolic.jl | 39 ++++++++++ src/manifolds/Sphere.jl | 42 +++++++++++ ...ymmetricPositiveDefiniteAffineInvariant.jl | 24 ++++++ src/utils.jl | 75 +++++++++++++++++++ test/manifolds/embedded_torus.jl | 22 +++--- test/manifolds/euclidean.jl | 6 +- test/manifolds/hyperbolic.jl | 33 ++++---- test/manifolds/sphere.jl | 13 +++- tutorials/_quarto.yml | 2 +- tutorials/exploring-curvature.qmd | 64 ++++++++++++++++ 18 files changed, 345 insertions(+), 42 deletions(-) create mode 100644 tutorials/exploring-curvature.qmd diff --git a/NEWS.md b/NEWS.md index 5484461eda..a00a23ad15 100644 --- a/NEWS.md +++ b/NEWS.md @@ -7,6 +7,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [0.9.15] unreleased +## Added + +* Implementation of `sectional_curvature`, `sectional_curvature_min` and `sectional_curvature_max` for several manifolds. +* `sectional_curvature_matrix` function and a tutorial on coordinate-free curvature. + ## Fixed * several typographical errors in the docs diff --git a/Project.toml b/Project.toml index 01f0dd47e7..3747090daa 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Manifolds" uuid = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e" authors = ["Seth Axen ", "Mateusz Baran ", "Ronny Bergmann ", "Antoine Levitt "] -version = "0.9.14" +version = "0.9.15" [deps] Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" @@ -51,7 +51,7 @@ HybridArrays = "0.4" Kronecker = "0.4, 0.5" LinearAlgebra = "1.6" ManifoldDiff = "0.3.7" -ManifoldsBase = "0.15.6" +ManifoldsBase = "0.15.8" Markdown = "1.6" MatrixEquations = "2.2" OrdinaryDiffEq = "6.31" diff --git a/docs/make.jl b/docs/make.jl index 74e80296e3..7c80a1f1f0 100755 --- a/docs/make.jl +++ b/docs/make.jl @@ -101,6 +101,7 @@ makedocs(; "work in charts" => "tutorials/working-in-charts.md", "perform Hand gesture analysis" => "tutorials/hand-gestures.md", "integrate on manifolds and handle probability densities" => "tutorials/integration.md", + "explore curvature without coordinates" => "tutorials/exploring-curvature.md", ], "Manifolds" => [ "Basic manifolds" => [ diff --git a/docs/src/features/utilities.md b/docs/src/features/utilities.md index d0fa97ca44..b2de5f384c 100644 --- a/docs/src/features/utilities.md +++ b/docs/src/features/utilities.md @@ -8,18 +8,12 @@ The following terms introduce a nicer notation for some operations, for example in ```` -## Fallback for the exponential map: Solving the corresponding ODE - -When additionally loading [`NLSolve.jl`](https://github.com/JuliaNLSolvers/NLsolve.jl) the following fallback for the exponential map is available. +## Public documentation -```@autodocs -Modules = [Manifolds] -Pages = ["nlsolve.jl"] -Order = [:type, :function] +```@docs +sectional_curvature_matrix ``` -## Public documentation - ### Specific exception types For some manifolds it is useful to keep an extra index, at which point on the manifold, the error occurred as well as to collect all errors that occurred on a manifold. This page contains the manifold-specific error messages this package introduces. diff --git a/docs/src/misc/notation.md b/docs/src/misc/notation.md index a87fd187a8..a1f0b393d5 100644 --- a/docs/src/misc/notation.md +++ b/docs/src/misc/notation.md @@ -49,6 +49,7 @@ Within the documented functions, the utf8 symbols are used whenever possible, as | ``\mathcal P_{t_1\gets t_0}^cX`` | parallel transport along the curve ``c``| ``\mathcal P^cX=\mathcal P_{1\gets 0}^cX`` | of the vector ``X`` from ``p=c(0)`` to ``c(1)`` | ``p`` | a point on ``\mathcal M`` | ``p_1, p_2, \ldots,q`` | for 3 points one might use ``x,y,z`` | | ``\operatorname{retr}_pX``| a retraction | | +| ``\kappa_p(X, Y)`` | sectional curvature | | | ``ξ`` | a set of tangent vectors | ``\{X_1,\ldots,X_n\}`` | | | ``J_{2n} \in ℝ^{2n×2n}`` | the [`SymplecticElement`](@ref) | | | | ``T_p \mathcal M`` | the tangent space at ``p`` | | | diff --git a/docs/src/references.bib b/docs/src/references.bib index 3f2cc96041..b0ef80ef9a 100644 --- a/docs/src/references.bib +++ b/docs/src/references.bib @@ -194,6 +194,16 @@ @article{LeBrignantPuechmorel:2019 # # C # ---------------------------------------------------------------------------------------- + +@book{CheegerEbin:2008, + address = {Providence, R.I}, + title = {Comparison {Theorems} in {Riemannian} {Geometry}}, + isbn = {978-0-8218-4417-5}, + publisher = {American Mathematical Society}, + author = {Cheeger, Jeffrey and Ebin, David G.}, + month = may, + year = {2008}, +} @book{Chikuse:2003, DOI = {10.1007/978-0-387-21540-2}, YEAR = {2003}, diff --git a/src/Manifolds.jl b/src/Manifolds.jl index c0d5b3e746..a4fde53b16 100644 --- a/src/Manifolds.jl +++ b/src/Manifolds.jl @@ -146,6 +146,9 @@ import ManifoldsBase: retract_softmax!, riemann_tensor, riemann_tensor!, + sectional_curvature, + sectional_curvature_max, + sectional_curvature_min, set_component!, submanifold, submanifold_component, @@ -906,6 +909,10 @@ export ×, riemannian_Hessian!, riemann_tensor, riemann_tensor!, + sectional_curvature, + sectional_curvature_matrix, + sectional_curvature_max, + sectional_curvature_min, set_component!, sharp, sharp!, diff --git a/src/manifolds/Euclidean.jl b/src/manifolds/Euclidean.jl index a0eacec744..861bdd6849 100644 --- a/src/manifolds/Euclidean.jl +++ b/src/manifolds/Euclidean.jl @@ -753,6 +753,33 @@ function riemann_tensor!(::Euclidean, Xresult, p, X, Y, Z) return fill!(Xresult, 0) end +@doc raw""" + sectional_curvature(::Euclidean, p, X, Y) + +Sectional curvature of [`Euclidean`](@ref) manifold `M` is 0. +""" +function sectional_curvature(::Euclidean, p, X, Y) + return 0.0 +end + +@doc raw""" + sectional_curvature_max(::Euclidean) + +Sectional curvature of [`Euclidean`](@ref) manifold `M` is 0. +""" +function sectional_curvature_max(::Euclidean) + return 0.0 +end + +@doc raw""" + sectional_curvature_min(M::Euclidean) + +Sectional curvature of [`Euclidean`](@ref) manifold `M` is 0. +""" +function sectional_curvature_min(::Euclidean) + return 0.0 +end + function Base.show(io::IO, M::Euclidean{N,𝔽}) where {N<:Tuple,𝔽} size = get_parameter(M.size) return print(io, "Euclidean($(join(size, ", ")); field=$(𝔽), parameter=:field)") diff --git a/src/manifolds/Hyperbolic.jl b/src/manifolds/Hyperbolic.jl index 8c1f55195e..e7d761ec2b 100644 --- a/src/manifolds/Hyperbolic.jl +++ b/src/manifolds/Hyperbolic.jl @@ -415,3 +415,42 @@ function riemann_tensor!(M::Hyperbolic, W, p, X, Y, Z) W .= inner(M, p, Z, X) .* Y .- inner(M, p, Z, Y) .* X return W end + +@doc raw""" + sectional_curvature(::Hyperbolic, p, X, Y) + +Sectional curvature of [`Hyperbolic`](@ref) `M` is -1 if dimension is > 1 and 0 otherwise. +""" +function sectional_curvature(M::Hyperbolic, p, X, Y) + if manifold_dimension(M) > 1 + return -1.0 + else + return 0.0 + end +end + +@doc raw""" + sectional_curvature_max(::Hyperbolic) + +Sectional curvature of [`Hyperbolic`](@ref) `M` is -1 if dimension is > 1 and 0 otherwise. +""" +function sectional_curvature_max(M::Hyperbolic) + if manifold_dimension(M) > 1 + return -1.0 + else + return 0.0 + end +end + +@doc raw""" + sectional_curvature_min(M::Hyperbolic) + +Sectional curvature of [`Hyperbolic`](@ref) `M` is -1 if dimension is > 1 and 0 otherwise. +""" +function sectional_curvature_min(M::Hyperbolic) + if manifold_dimension(M) > 1 + return -1.0 + else + return 0.0 + end +end diff --git a/src/manifolds/Sphere.jl b/src/manifolds/Sphere.jl index 8ed626caaf..0c76c085fd 100644 --- a/src/manifolds/Sphere.jl +++ b/src/manifolds/Sphere.jl @@ -571,6 +571,48 @@ function riemann_tensor!(M::AbstractSphere{ℝ}, Xresult, p, X, Y, Z) return Xresult end +@doc raw""" + sectional_curvature(::AbstractSphere, p, X, Y) + +Sectional curvature of [`AbstractSphere`](@ref) `M` is 1 if dimension is greater than 1 +and 0 otherwise. +""" +function sectional_curvature(M::AbstractSphere, p, X, Y) + if manifold_dimension(M) > 1 + return 1.0 + else + return 0.0 + end +end + +@doc raw""" + sectional_curvature_max(::AbstractSphere) + +Sectional curvature of [`AbstractSphere`](@ref) `M` is 1 if dimension is greater than 1 +and 0 otherwise. +""" +function sectional_curvature_max(M::AbstractSphere) + if manifold_dimension(M) > 1 + return 1.0 + else + return 0.0 + end +end + +@doc raw""" + sectional_curvature_min(M::AbstractSphere) + +Sectional curvature of [`AbstractSphere`](@ref) `M` is 1 if dimension is greater than 1 +and 0 otherwise. +""" +function sectional_curvature_min(M::AbstractSphere) + if manifold_dimension(M) > 1 + return 1.0 + else + return 0.0 + end +end + @doc raw""" volume_density(M::AbstractSphere{ℝ}, p, X) diff --git a/src/manifolds/SymmetricPositiveDefiniteAffineInvariant.jl b/src/manifolds/SymmetricPositiveDefiniteAffineInvariant.jl index 1c171e93b2..0500415530 100644 --- a/src/manifolds/SymmetricPositiveDefiniteAffineInvariant.jl +++ b/src/manifolds/SymmetricPositiveDefiniteAffineInvariant.jl @@ -454,6 +454,30 @@ function riemann_tensor!(::SymmetricPositiveDefinite, Xresult, p, X, Y, Z) return Xresult end +""" + sectional_curvature_min(M::SymmetricPositiveDefinite) + +Return minimum sectional curvature of [`SymmetricPositiveDefinite`](@ref) manifold, +that is 0 for SPD(1) and SPD(2) and -0.25 otherwise. +""" +function sectional_curvature_min(M::SymmetricPositiveDefinite) + if manifold_dimension(M) < 2 + return 0.0 + else + return -0.25 + end +end + +""" + sectional_curvature_max(M::SymmetricPositiveDefinite) + +Return minimum sectional curvature of [`SymmetricPositiveDefinite`](@ref) manifold, +that is 0. +""" +function sectional_curvature_max(::SymmetricPositiveDefinite) + return 0.0 +end + """ volume_density(::SymmetricPositiveDefinite, p, X) diff --git a/src/utils.jl b/src/utils.jl index 4dcf8757e4..c5f5873729 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -277,3 +277,78 @@ _eps_safe(::Type{T}) where {T<:Real} = eps(T) _eps_safe(::Type{T}) where {T<:Number} = eps(real(T)) max_eps(xs...) = maximum(_eps_safe ∘ eltype, xs) + +""" + sectional_curvature_matrix(M::AbstractManifold, p, B::AbstractBasis) + +Compute the matrix of sectional curvatures of manifold `M` at point `p`. +Entry `(i, j)` corresponds to sectional curvature of the surface spanned by vectors +`i` and `j` from basis `B`. +""" +function sectional_curvature_matrix(M::AbstractManifold, p, B::AbstractBasis) + V = get_vectors(M, p, get_basis(M, p, B)) + N = length(V) + result = zeros(N, N) + for (i, e_i) in enumerate(V) + for (j, e_j) in enumerate(V) + if i < j + result[i, j] = sectional_curvature(M, p, e_i, e_j) + result[j, i] = result[i, j] + end + end + end + return result +end + +@doc raw""" + sectional_curvature_bdp(M::AbstractManifold, p, X, Y; r::Real=1e-3, N::Int=10000) + +Approximate sectional curvature of manifold `M` in the plane spanned by vectors `X` and `Y` +from tangent space at `p` using a circle on `M` of radius `r` divided into `N` segments. + +The approximation is derived from the [Bertrand–Diguet–Puiseux theorem](https://en.wikipedia.org/wiki/Bertrand%E2%80%93Diguet%E2%80%93Puiseux_theorem) +which states that +````math +\kappa_p(X, Y) = \lim_{r \to 0^+} 3\frac{2\pi r-C(r)}{\pi r^3}, +```` +where ``C(r)`` is the circumference of the circle of radius ``r`` around `p` in submanifold +of `M` spanned by `X` and `Y`. The circumference calculation method has a tendency to +return curvature values larger than the exact ones. +""" +function sectional_curvature_bdp(M::AbstractManifold, p, X, Y; r::Real=1e-3, N::Int=10000) + circumference = 0.0 + p_i = similar(p) + p_ip1 = similar(p) + for i in 1:N + θ_i = 2π * (i-1) / N + θ_ip1 = 2π * (i) / N + exp!(M, p_i, p, r .* (sin(θ_i) .* X .+ cos(θ_i) .* Y)) + exp!(M, p_ip1, p, r .* (sin(θ_ip1) .* X .+ cos(θ_ip1) .* Y)) + + circumference += distance(M, p_i, p_ip1) + end + return 3 * (2π * r - circumference) / (π * r^3) +end + +""" + sectional_curvature_matrix_bdp(M::AbstractManifold, p, B::AbstractBasis; r::Real=1e-3, N::Int=10000) + +Estimate the matrix of sectional curvatures of manifold `M` at point `p` using +`sectional_curvature_bdp`. Entry `(i, j)`` corresponds to sectional curvature of the +surface spanned by vectors `i` and `j` from basis `B`. +""" +function sectional_curvature_matrix_bdp(M::AbstractManifold, p, B::AbstractBasis; r::Real=1e-3, N_pts::Int=10000) + V = get_vectors(M, p, get_basis(M, p, B)) + N = length(V) + result = zeros(N, N) + for (i, e_i) in enumerate(V) + for (j, e_j) in enumerate(V) + if i < j + result[i, j] = sectional_curvature_bdp(M, p, e_i, e_j; r=r, N=N_pts) + result[j, i] = result[i, j] + end + end + end + return result +end + diff --git a/test/manifolds/embedded_torus.jl b/test/manifolds/embedded_torus.jl index b8720d10d5..1916351dc3 100644 --- a/test/manifolds/embedded_torus.jl +++ b/test/manifolds/embedded_torus.jl @@ -89,17 +89,17 @@ using BoundaryValueDiffEq Manifolds.transition_map_diff(M, A, i_p0x, [0.0, 0.0], X_p0x, (-1.0, -0.3)) a2 = [-0.5, 0.3] - sol_log = Manifolds.solve_chart_log_bvp(M, p0x, a2, A, (0, 0)) - @test sol_log(0.0)[1:2] ≈ p0x - @test sol_log(1.0)[1:2] ≈ a2 atol = 1e-7 - # a test randomly failed here on Julia 1.6 once for no clear reason? - # so I bumped tolerance considerably - bvp_atol = VERSION < v"1.7" ? 2e-3 : 1e-15 - @test isapprox( - norm(M, A, (0, 0), p0x, sol_log(0.0)[3:4]), - Manifolds.estimate_distance_from_bvp(M, p0x, a2, A, (0, 0)); - atol=bvp_atol, - ) + # sol_log = Manifolds.solve_chart_log_bvp(M, p0x, a2, A, (0, 0)) + # @test sol_log(0.0)[1:2] ≈ p0x + # @test sol_log(1.0)[1:2] ≈ a2 atol = 1e-7 + # # a test randomly failed here on Julia 1.6 once for no clear reason? + # # so I bumped tolerance considerably + # bvp_atol = VERSION < v"1.7" ? 2e-3 : 1e-15 + # @test isapprox( + # norm(M, A, (0, 0), p0x, sol_log(0.0)[3:4]), + # Manifolds.estimate_distance_from_bvp(M, p0x, a2, A, (0, 0)); + # atol=bvp_atol, + # ) @test Manifolds.IntegratorTerminatorNearChartBoundary().check_chart_switch_kwargs === NamedTuple() diff --git a/test/manifolds/euclidean.jl b/test/manifolds/euclidean.jl index 64e856177a..0099f030ec 100644 --- a/test/manifolds/euclidean.jl +++ b/test/manifolds/euclidean.jl @@ -296,8 +296,12 @@ using FiniteDifferences @test size(RT) == (2, 2, 2, 2) @test norm(RT) ≈ 0.0 atol = 1e-16 - @test riemann_tensor(M, p, [1, 2], [1, 3], [1, 4]) == [0, 0] @test !Manifolds.check_chart_switch(M, A, i, p) + + @test riemann_tensor(M, p, [1, 2], [1, 3], [1, 4]) == [0, 0] + @test sectional_curvature(M, p, [1.0, 0.0], [0.0, 1.0]) == 0.0 + @test sectional_curvature_max(M) == 0.0 + @test sectional_curvature_min(M) == 0.0 end @testset "Induced Basis and local metric for EuclideanMetric" begin struct DefaultManifold <: AbstractManifold{ℝ} end diff --git a/test/manifolds/hyperbolic.jl b/test/manifolds/hyperbolic.jl index f4970ca3d1..07e18b4d67 100644 --- a/test/manifolds/hyperbolic.jl +++ b/test/manifolds/hyperbolic.jl @@ -347,23 +347,22 @@ include("../header.jl") end @testset "other metric" begin M = Hyperbolic(2) - @test riemann_tensor( - M, - [0.0, 0.0, -1.0], - [-1.0, 0.0, 0.0], - [0.0, -1.0, 0.0], - [-0.5, -0.5, 0.0], - ) == [0.5, -0.5, 0.0] - X = [0.0, 0.0, 0.0] - @test riemann_tensor!( - M, - X, - [0.0, 0.0, -1.0], - [-1.0, 0.0, 0.0], - [0.0, -1.0, 0.0], - [-0.5, -0.5, 0.0], - ) === X - @test X == [0.5, -0.5, 0.0] + p = [0.0, 0.0, -1.0] + X = [-1.0, 0.0, 0.0] + Y = [0.0, -1.0, 0.0] + Z = [-0.5, -0.5, 0.0] + @test riemann_tensor(M, p, X, Y, Z) == [0.5, -0.5, 0.0] + Xout = [0.0, 0.0, 0.0] + @test riemann_tensor!(M, Xout, p, X, Y, Z) === Xout + @test Xout == [0.5, -0.5, 0.0] + + @test sectional_curvature(M, p, X, Y) == -1.0 + @test sectional_curvature_max(M) == -1.0 + @test sectional_curvature_min(M) == -1.0 + M1 = Hyperbolic(1) + @test sectional_curvature(M1, p, X, Y) == 0.0 + @test sectional_curvature_max(M1) == 0.0 + @test sectional_curvature_min(M1) == 0.0 end @testset "ManifoldDiff" begin # ManifoldDiff diff --git a/test/manifolds/sphere.jl b/test/manifolds/sphere.jl index 9031742fb2..e81c340079 100644 --- a/test/manifolds/sphere.jl +++ b/test/manifolds/sphere.jl @@ -255,7 +255,18 @@ using ManifoldsBase: TFVector @testset "other metric" begin M = Sphere(2) - @test riemann_tensor(M, [0, 1, 0], [1, 0, 0], [1, 0, 4], [0, 0, 1]) == [4, 0, 0] + p = [0, 1, 0] + X = [1, 0, 0] + Y = [1, 0, 4] + Z = [0, 0, 1] + @test riemann_tensor(M, p, X, Y, Z) == [4, 0, 0] + @test sectional_curvature(M, p, X, Y) == 1.0 + @test sectional_curvature_max(M) == 1.0 + @test sectional_curvature_min(M) == 1.0 + M1 = Sphere(1) + @test sectional_curvature(M1, p, X, Y) == 0.0 + @test sectional_curvature_max(M1) == 0.0 + @test sectional_curvature_min(M1) == 0.0 end @testset "ManifoldDiff" begin diff --git a/tutorials/_quarto.yml b/tutorials/_quarto.yml index d046c3b30f..5befbaf19d 100644 --- a/tutorials/_quarto.yml +++ b/tutorials/_quarto.yml @@ -24,4 +24,4 @@ format: variant: -raw_html+tex_math_dollars wrap: preserve -jupyter: julia-1.9 \ No newline at end of file +jupyter: julia-1.10 diff --git a/tutorials/exploring-curvature.qmd b/tutorials/exploring-curvature.qmd new file mode 100644 index 0000000000..6e65d4424f --- /dev/null +++ b/tutorials/exploring-curvature.qmd @@ -0,0 +1,64 @@ +--- +title: Exploring curvature without coordinates +--- + +```{julia} +#| echo: false +#| code-fold: true +#| output: false +using Pkg; +cd(@__DIR__) +Pkg.activate("."); # for reproducibility use the local tutorial environment. +Pkg.develop(PackageSpec(; path=(@__DIR__) * "/../")) +using Markdown +``` + +This part of documentation covers exploration of curvature of manifolds $\mathcal{M}$. +There are multiple ways to describe curvature: Christoffel symbols, Riemann tensor, Ricci tensor, sectional curvature, and many other. +They are usually considered only in coordinates but there is a way to demonstrate curvature in coordinate-free way. + +## Sectional curvature matrix + +Curvature of a manifold can be explored using the [`sectional_curvature_matrix`](@ref) function. +Note that Riemann tensor and sectional curvature are equivalently full specifications of curvature in a manifold, see [CheegerEbin:2008](@cite), Eq. (1.12). + +Let's take the [`SymmetricPositiveDefinite`](@ref) manifold as our first example. +It has nonpositive sectional curvature: + +```{julia} +using Manifolds +M = SymmetricPositiveDefinite(3) +p = rand(M) +cm = sectional_curvature_matrix(M, p, DefaultOrthonormalBasis()) +``` + +We can verify that the curvature is consistent with an approximation based on the Bertrand–Diguet–Puiseux theorem, which relies only on an ONB, exponential map and distance calculation: +```{julia} +cm_bdp = Manifolds.sectional_curvature_matrix_bdp(M, p, DefaultOrthonormalBasis(); r=1e-3, N_pts=100000) +println(norm(cm - cm_bdp)) +``` + +This approximation converges quite slowly with `N_pts` and is prone to numerical errors at low values of `r` and large values of `N_pts`. + +You can also take the vectors from the basis and see what kind of planes they correspond to. +It may be easier to see for the identity matrix as the base point. +```{julia} +p = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0] +V = get_vectors(M, p, get_basis(M, p, DefaultOrthonormalBasis())) +cm = sectional_curvature_matrix(M, p, DefaultOrthonormalBasis()) +for X in V + println(exp(M, p, X)) +end +``` + +The flat planes correspond to directions where the matrix changes independently. +In other cases sectional curvature indicates hyperbolic characteristic of a submanifold. + + + +## Literature + +```@bibliography +Pages = ["exploring-curvature.md"] +Canonical=false +``` \ No newline at end of file From 8a1e092b2603695531d7cb24ec23c7ac35513ae0 Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Sat, 16 Mar 2024 09:53:27 +0100 Subject: [PATCH 02/14] formatting --- src/utils.jl | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/src/utils.jl b/src/utils.jl index c5f5873729..d09d05b20d 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -320,7 +320,7 @@ function sectional_curvature_bdp(M::AbstractManifold, p, X, Y; r::Real=1e-3, N:: p_i = similar(p) p_ip1 = similar(p) for i in 1:N - θ_i = 2π * (i-1) / N + θ_i = 2π * (i - 1) / N θ_ip1 = 2π * (i) / N exp!(M, p_i, p, r .* (sin(θ_i) .* X .+ cos(θ_i) .* Y)) exp!(M, p_ip1, p, r .* (sin(θ_ip1) .* X .+ cos(θ_ip1) .* Y)) @@ -337,7 +337,13 @@ Estimate the matrix of sectional curvatures of manifold `M` at point `p` using `sectional_curvature_bdp`. Entry `(i, j)`` corresponds to sectional curvature of the surface spanned by vectors `i` and `j` from basis `B`. """ -function sectional_curvature_matrix_bdp(M::AbstractManifold, p, B::AbstractBasis; r::Real=1e-3, N_pts::Int=10000) +function sectional_curvature_matrix_bdp( + M::AbstractManifold, + p, + B::AbstractBasis; + r::Real=1e-3, + N_pts::Int=10000, +) V = get_vectors(M, p, get_basis(M, p, B)) N = length(V) result = zeros(N, N) @@ -351,4 +357,3 @@ function sectional_curvature_matrix_bdp(M::AbstractManifold, p, B::AbstractBasis end return result end - From c5de3ab17b2526fa4186d19e300c2ae6d09f0524 Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Sat, 16 Mar 2024 10:19:47 +0100 Subject: [PATCH 03/14] update CI settings --- .github/workflows/documenter.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/documenter.yml b/.github/workflows/documenter.yml index 968830a64c..aee47acf41 100644 --- a/.github/workflows/documenter.yml +++ b/.github/workflows/documenter.yml @@ -14,10 +14,10 @@ jobs: - uses: actions/checkout@v4 - uses: quarto-dev/quarto-actions/setup@v2 with: - version: 1.3.361 + version: 1.4.551 - uses: julia-actions/setup-julia@latest with: - version: 1.9 + version: 1.10 - name: Julia Cache uses: julia-actions/cache@v1 - name: Cache Quarto From b961c2c3983d5263eb17ef0755f43c83e81ddf5b Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Sat, 16 Mar 2024 10:23:02 +0100 Subject: [PATCH 04/14] fix again --- .github/workflows/documenter.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/documenter.yml b/.github/workflows/documenter.yml index aee47acf41..3c9d5d4372 100644 --- a/.github/workflows/documenter.yml +++ b/.github/workflows/documenter.yml @@ -14,10 +14,10 @@ jobs: - uses: actions/checkout@v4 - uses: quarto-dev/quarto-actions/setup@v2 with: - version: 1.4.551 + version: "1.4.551" - uses: julia-actions/setup-julia@latest with: - version: 1.10 + version: "1.10" - name: Julia Cache uses: julia-actions/cache@v1 - name: Cache Quarto From 72b9aafe6b2eca47f37e4fa2da4cc4bc9e0433e8 Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Wed, 20 Mar 2024 09:05:54 +0100 Subject: [PATCH 05/14] address review --- src/utils.jl | 19 +++++++++++++------ tutorials/exploring-curvature.qmd | 2 +- 2 files changed, 14 insertions(+), 7 deletions(-) diff --git a/src/utils.jl b/src/utils.jl index d09d05b20d..5c31fd2038 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -301,7 +301,7 @@ function sectional_curvature_matrix(M::AbstractManifold, p, B::AbstractBasis) end @doc raw""" - sectional_curvature_bdp(M::AbstractManifold, p, X, Y; r::Real=1e-3, N::Int=10000) + estimated_sectional_curvature(M::AbstractManifold, p, X, Y; r::Real=1e-3, N::Int=10000) Approximate sectional curvature of manifold `M` in the plane spanned by vectors `X` and `Y` from tangent space at `p` using a circle on `M` of radius `r` divided into `N` segments. @@ -315,7 +315,14 @@ where ``C(r)`` is the circumference of the circle of radius ``r`` around `p` in of `M` spanned by `X` and `Y`. The circumference calculation method has a tendency to return curvature values larger than the exact ones. """ -function sectional_curvature_bdp(M::AbstractManifold, p, X, Y; r::Real=1e-3, N::Int=10000) +function estimated_sectional_curvature( + M::AbstractManifold, + p, + X, + Y; + r::Real=1e-3, + N::Int=10000, +) circumference = 0.0 p_i = similar(p) p_ip1 = similar(p) @@ -331,13 +338,13 @@ function sectional_curvature_bdp(M::AbstractManifold, p, X, Y; r::Real=1e-3, N:: end """ - sectional_curvature_matrix_bdp(M::AbstractManifold, p, B::AbstractBasis; r::Real=1e-3, N::Int=10000) + estimated_sectional_curvature_matrix(M::AbstractManifold, p, B::AbstractBasis; r::Real=1e-3, N::Int=10000) Estimate the matrix of sectional curvatures of manifold `M` at point `p` using -`sectional_curvature_bdp`. Entry `(i, j)`` corresponds to sectional curvature of the +`estimated_sectional_curvature`. Entry `(i, j)`` corresponds to sectional curvature of the surface spanned by vectors `i` and `j` from basis `B`. """ -function sectional_curvature_matrix_bdp( +function estimated_sectional_curvature_matrix( M::AbstractManifold, p, B::AbstractBasis; @@ -350,7 +357,7 @@ function sectional_curvature_matrix_bdp( for (i, e_i) in enumerate(V) for (j, e_j) in enumerate(V) if i < j - result[i, j] = sectional_curvature_bdp(M, p, e_i, e_j; r=r, N=N_pts) + result[i, j] = estimated_sectional_curvature(M, p, e_i, e_j; r=r, N=N_pts) result[j, i] = result[i, j] end end diff --git a/tutorials/exploring-curvature.qmd b/tutorials/exploring-curvature.qmd index 6e65d4424f..a800385275 100644 --- a/tutorials/exploring-curvature.qmd +++ b/tutorials/exploring-curvature.qmd @@ -34,7 +34,7 @@ cm = sectional_curvature_matrix(M, p, DefaultOrthonormalBasis()) We can verify that the curvature is consistent with an approximation based on the Bertrand–Diguet–Puiseux theorem, which relies only on an ONB, exponential map and distance calculation: ```{julia} -cm_bdp = Manifolds.sectional_curvature_matrix_bdp(M, p, DefaultOrthonormalBasis(); r=1e-3, N_pts=100000) +cm_bdp = Manifolds.estimated_sectional_curvature_matrix(M, p, DefaultOrthonormalBasis(); r=1e-3, N_pts=100000) println(norm(cm - cm_bdp)) ``` From c282da84052d0b0dbb7d44ddf6e3e074e5570dd9 Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Wed, 20 Mar 2024 10:02:20 +0100 Subject: [PATCH 06/14] Upper-bounding NonlinearSolve --- Project.toml | 2 ++ test/manifolds/embedded_torus.jl | 22 +++++++++++----------- 2 files changed, 13 insertions(+), 11 deletions(-) diff --git a/Project.toml b/Project.toml index 3747090daa..0a0886f218 100644 --- a/Project.toml +++ b/Project.toml @@ -54,6 +54,7 @@ ManifoldDiff = "0.3.7" ManifoldsBase = "0.15.8" Markdown = "1.6" MatrixEquations = "2.2" +NonlinearSolve = "3 - 3.7.3" OrdinaryDiffEq = "6.31" Plots = "1" Quaternions = "0.5, 0.6, 0.7" @@ -78,6 +79,7 @@ Gtk = "4c0ca9eb-093a-5379-98c5-f87ac0bbbf44" ImageIO = "82e4d734-157c-48bb-816b-45c225c6df19" ImageMagick = "6218d12a-5da1-5696-b52f-db25d2ecc6d1" NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" +NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" PythonPlot = "274fc56d-3b97-40fa-a1cd-1b4a50311bf9" diff --git a/test/manifolds/embedded_torus.jl b/test/manifolds/embedded_torus.jl index 1916351dc3..b8720d10d5 100644 --- a/test/manifolds/embedded_torus.jl +++ b/test/manifolds/embedded_torus.jl @@ -89,17 +89,17 @@ using BoundaryValueDiffEq Manifolds.transition_map_diff(M, A, i_p0x, [0.0, 0.0], X_p0x, (-1.0, -0.3)) a2 = [-0.5, 0.3] - # sol_log = Manifolds.solve_chart_log_bvp(M, p0x, a2, A, (0, 0)) - # @test sol_log(0.0)[1:2] ≈ p0x - # @test sol_log(1.0)[1:2] ≈ a2 atol = 1e-7 - # # a test randomly failed here on Julia 1.6 once for no clear reason? - # # so I bumped tolerance considerably - # bvp_atol = VERSION < v"1.7" ? 2e-3 : 1e-15 - # @test isapprox( - # norm(M, A, (0, 0), p0x, sol_log(0.0)[3:4]), - # Manifolds.estimate_distance_from_bvp(M, p0x, a2, A, (0, 0)); - # atol=bvp_atol, - # ) + sol_log = Manifolds.solve_chart_log_bvp(M, p0x, a2, A, (0, 0)) + @test sol_log(0.0)[1:2] ≈ p0x + @test sol_log(1.0)[1:2] ≈ a2 atol = 1e-7 + # a test randomly failed here on Julia 1.6 once for no clear reason? + # so I bumped tolerance considerably + bvp_atol = VERSION < v"1.7" ? 2e-3 : 1e-15 + @test isapprox( + norm(M, A, (0, 0), p0x, sol_log(0.0)[3:4]), + Manifolds.estimate_distance_from_bvp(M, p0x, a2, A, (0, 0)); + atol=bvp_atol, + ) @test Manifolds.IntegratorTerminatorNearChartBoundary().check_chart_switch_kwargs === NamedTuple() From f64fc579e0ec2f390683037077d919d40bb753d4 Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Wed, 20 Mar 2024 10:36:38 +0100 Subject: [PATCH 07/14] Add NonlinearSolve as a direct dep --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 0a0886f218..c68d527dd2 100644 --- a/Project.toml +++ b/Project.toml @@ -14,6 +14,7 @@ ManifoldDiff = "af67fdf4-a580-4b9f-bbec-742ef357defd" ManifoldsBase = "3362f125-f0bb-47a3-aa74-596ffd7ef2fb" Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a" MatrixEquations = "99c1a7ee-ab34-5fd5-8076-27c950a045f4" +NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" Quaternions = "94ee1d12-ae83-5a48-8b1c-48b8ff168ae0" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" @@ -79,7 +80,6 @@ Gtk = "4c0ca9eb-093a-5379-98c5-f87ac0bbbf44" ImageIO = "82e4d734-157c-48bb-816b-45c225c6df19" ImageMagick = "6218d12a-5da1-5696-b52f-db25d2ecc6d1" NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" -NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" PythonPlot = "274fc56d-3b97-40fa-a1cd-1b4a50311bf9" From 9c1a8f6c0ce69d61480fe767f7055784bbbd87b4 Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Wed, 20 Mar 2024 10:42:46 +0100 Subject: [PATCH 08/14] tweaking --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index c68d527dd2..4ba1e082c6 100644 --- a/Project.toml +++ b/Project.toml @@ -55,7 +55,7 @@ ManifoldDiff = "0.3.7" ManifoldsBase = "0.15.8" Markdown = "1.6" MatrixEquations = "2.2" -NonlinearSolve = "3 - 3.7.3" +NonlinearSolve = "1 - 3.7.3" OrdinaryDiffEq = "6.31" Plots = "1" Quaternions = "0.5, 0.6, 0.7" From 1f10463976732bfcf85c531469f06d2f5517cc70 Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Thu, 21 Mar 2024 11:43:01 +0100 Subject: [PATCH 09/14] more curvature --- docs/src/references.bib | 12 ++++++++++++ src/manifolds/Rotations.jl | 29 +++++++++++++++++++++++++++++ test/manifolds/rotations.jl | 8 ++++++++ tutorials/exploring-curvature.qmd | 31 +++++++++++++++++++++++++++++++ 4 files changed, 80 insertions(+) diff --git a/docs/src/references.bib b/docs/src/references.bib index b0ef80ef9a..e377fd0737 100644 --- a/docs/src/references.bib +++ b/docs/src/references.bib @@ -351,6 +351,18 @@ @article{GaoSonAbsilStykel:2021 TITLE = {Riemannian Optimization on the Symplectic Stiefel Manifold}, JOURNAL = {SIAM Journal on Optimization} } +@article{Ge:2014, + title = {{DDVV}-type inequality for skew-symmetric matrices and {Simons}-type inequality for {Riemannian} submersions}, + volume = {251}, + issn = {0001-8708}, + url = {https://www.sciencedirect.com/science/article/pii/S000187081300385X}, + doi = {10.1016/j.aim.2013.10.010}, + journal = {Advances in Mathematics}, + author = {Ge, Jianquan}, + month = jan, + year = {2014}, + pages = {62--86}, +} @inproceedings{Giles:2008, address = {Berlin, Heidelberg}, series = {Lecture {Notes} in {Computational} {Science} and {Engineering}}, diff --git a/src/manifolds/Rotations.jl b/src/manifolds/Rotations.jl index dc66b52b7b..6f756bf9bc 100644 --- a/src/manifolds/Rotations.jl +++ b/src/manifolds/Rotations.jl @@ -423,6 +423,35 @@ function riemannian_Hessian!(M::Rotations, Y, p, G, H, X) return Y end +@doc raw""" + sectional_curvature_max(::Rotations) + +Sectional curvature of [`Rotations`](@ref) `M` is equal to 0 for `Rotations(1)` and +`Rotations(2)`, less than or equal to 1/8 for `Rotations(3)` and less than or equal to 1/4 +for higher-dimensional rotations manifolds. + +For reference, see [Ge:2014](@cite), Lemma 2.5 and [CheegerEbin:2008](@cite), Corollary 3.19. +""" +function sectional_curvature_max(M::Rotations) + N = manifold_dimension(M) + if N <= 1 + return 0.0 + elseif N == 3 + return 1 / 8 + else + return 1 / 4 + end +end + +@doc raw""" + sectional_curvature_min(M::Rotations) + +Sectional curvature of [`Rotations`](@ref) `M` is greater than or equal to 0. +""" +function sectional_curvature_min(::Rotations) + return 0.0 +end + Distributions.support(d::NormalRotationDistribution) = MPointSupport(d.manifold) @doc raw""" diff --git a/test/manifolds/rotations.jl b/test/manifolds/rotations.jl index dfd4aee46e..4b578a4687 100644 --- a/test/manifolds/rotations.jl +++ b/test/manifolds/rotations.jl @@ -275,6 +275,14 @@ include("../header.jl") ] end + @testset "sectional curvature" begin + @test sectional_curvature_min(Rotations(3)) == 0.0 + @test sectional_curvature_max(Rotations(1)) == 0.0 + @test sectional_curvature_max(Rotations(2)) == 0.0 + @test sectional_curvature_max(Rotations(3)) == 1 / 8 + @test sectional_curvature_max(Rotations(4)) == 1 / 4 + end + @testset "field parameter" begin M = Rotations(2; parameter=:field) @test is_flat(M) diff --git a/tutorials/exploring-curvature.qmd b/tutorials/exploring-curvature.qmd index a800385275..8d13e4f03c 100644 --- a/tutorials/exploring-curvature.qmd +++ b/tutorials/exploring-curvature.qmd @@ -27,6 +27,7 @@ It has nonpositive sectional curvature: ```{julia} using Manifolds +using LinearAlgebra M = SymmetricPositiveDefinite(3) p = rand(M) cm = sectional_curvature_matrix(M, p, DefaultOrthonormalBasis()) @@ -54,7 +55,37 @@ end The flat planes correspond to directions where the matrix changes independently. In other cases sectional curvature indicates hyperbolic characteristic of a submanifold. +Sectional curvature can be either larger or smaller than entries in the matrix on other planes. +Consider for example the manifold of rotation matrices in four dimensions, and a function that computes plane of maximum curvature using random search. +```{julia} +function max_curvature(M::AbstractManifold, p) + mc = -Inf + X = zero_vector(M, p) + Y = zero_vector(M, p) + for _ in 1:10000 + X_c = rand(M; vector_at=p) + Y_c = rand(M; vector_at=p) + sc = sectional_curvature(M, p, X_c, Y_c) + if sc > mc + mc = sc + X .= X_c + Y .= Y_c + end + end + return mc, X, Y +end + +M = Rotations(4) +p = Matrix(I(4) * 1.0) +println(sectional_curvature_matrix(M, p, DefaultOrthonormalBasis())) +mc, X, Y = max_curvature(M, p) +println(mc) +println(X) +println(Y) +``` + +In the planes corresponding to orthonormal basis, the maximum sectional curvature is 0.125 but the true upper bound is 0.25. ## Literature From 47de435d96cad637b84b7ba5887e1726ef29e854 Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Thu, 21 Mar 2024 12:57:27 +0100 Subject: [PATCH 10/14] add atol to failing test --- test/manifolds/multinomial_spd.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/manifolds/multinomial_spd.jl b/test/manifolds/multinomial_spd.jl index 70a6027a22..ea37d06023 100644 --- a/test/manifolds/multinomial_spd.jl +++ b/test/manifolds/multinomial_spd.jl @@ -34,6 +34,6 @@ include("../header.jl") @testset "Random" begin q = zeros(3, 3) M = MultinomialSymmetricPositiveDefinite(3) - @test is_point(M, rand!(MersenneTwister(), M, q)) + @test is_point(M, rand!(MersenneTwister(), M, q); atol=1e-15) end end From bbc6d11edf64b3d69687cfb8615ea61f52a5707d Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Thu, 21 Mar 2024 14:58:10 +0100 Subject: [PATCH 11/14] bump tolerance --- test/groups/power_group.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/groups/power_group.jl b/test/groups/power_group.jl index d2aa761b94..ede5304e18 100644 --- a/test/groups/power_group.jl +++ b/test/groups/power_group.jl @@ -30,7 +30,7 @@ include("group_utils.jl") pts, X_pts, X_pts; - atol=2e-9, + atol=2e-8, test_diff=true, test_log_from_identity=true, test_exp_from_identity=true, From b82c465e0fcd596faee55113047255e06555276e Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Sat, 23 Mar 2024 14:28:55 +0100 Subject: [PATCH 12/14] upper bound for PythonCall too --- Project.toml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Project.toml b/Project.toml index 4ba1e082c6..f9d8718cd4 100644 --- a/Project.toml +++ b/Project.toml @@ -16,6 +16,7 @@ Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a" MatrixEquations = "99c1a7ee-ab34-5fd5-8076-27c950a045f4" NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" Quaternions = "94ee1d12-ae83-5a48-8b1c-48b8ff168ae0" +PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" Requires = "ae029012-a4dd-5104-9daa-d747884805df" @@ -58,6 +59,7 @@ MatrixEquations = "2.2" NonlinearSolve = "1 - 3.7.3" OrdinaryDiffEq = "6.31" Plots = "1" +PythonCall = "0.9 - 0.9.15" Quaternions = "0.5, 0.6, 0.7" Random = "1.6" RecipesBase = "1.1" From 613e2831fbc77db25df463664d131bc6a1fb64ec Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Sat, 23 Mar 2024 16:20:51 +0100 Subject: [PATCH 13/14] change default vector transport method for some unitary matrix manifolds --- NEWS.md | 8 ++++++-- src/manifolds/GeneralUnitaryMatrices.jl | 2 ++ src/manifolds/Rotations.jl | 2 ++ test/manifolds/rotations.jl | 1 + test/manifolds/unitary_matrices.jl | 1 + 5 files changed, 12 insertions(+), 2 deletions(-) diff --git a/NEWS.md b/NEWS.md index a00a23ad15..65f13b660f 100644 --- a/NEWS.md +++ b/NEWS.md @@ -7,12 +7,16 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [0.9.15] unreleased -## Added +### Added * Implementation of `sectional_curvature`, `sectional_curvature_min` and `sectional_curvature_max` for several manifolds. * `sectional_curvature_matrix` function and a tutorial on coordinate-free curvature. -## Fixed +### Changed + +* `default_vector_transport_method` for `GeneralUnitaryMatrices` other than `Rotations` was changed to `ProjectionTransport`. + +### Fixed * several typographical errors in the docs * unifies to use two backticks ``` `` ``` for math instead of ` $ ` further in the docs diff --git a/src/manifolds/GeneralUnitaryMatrices.jl b/src/manifolds/GeneralUnitaryMatrices.jl index 7adcc0c002..db8b489bf2 100644 --- a/src/manifolds/GeneralUnitaryMatrices.jl +++ b/src/manifolds/GeneralUnitaryMatrices.jl @@ -182,6 +182,8 @@ function default_approximation_method(::GeneralUnitaryMatrices{<:Any,ℝ}, ::typ return GeodesicInterpolationWithinRadius(π / 2 / √2) end +default_vector_transport_method(::GeneralUnitaryMatrices) = ProjectionTransport() + embed(::GeneralUnitaryMatrices, p) = p @doc raw""" diff --git a/src/manifolds/Rotations.jl b/src/manifolds/Rotations.jl index 6f756bf9bc..ef3aee2cad 100644 --- a/src/manifolds/Rotations.jl +++ b/src/manifolds/Rotations.jl @@ -62,6 +62,8 @@ function angles_4d_skew_sym_matrix(A) return sqrt(halfb + sqrtdisc), sqrt(halfb - sqrtdisc) end +default_vector_transport_method(::Rotations) = ParallelTransport() + # from https://github.com/JuliaManifolds/Manifolds.jl/issues/453#issuecomment-1046057557 function _get_tridiagonal_elements(trian) N = size(trian, 1) diff --git a/test/manifolds/rotations.jl b/test/manifolds/rotations.jl index 4b578a4687..a36b646597 100644 --- a/test/manifolds/rotations.jl +++ b/test/manifolds/rotations.jl @@ -17,6 +17,7 @@ include("../header.jl") TEST_FLOAT32 && push!(types, Matrix{Float32}) TEST_STATIC_SIZED && push!(types, MMatrix{2,2,Float64,4}) retraction_methods = [PolarRetraction(), QRRetraction()] + @test default_vector_transport_method(M) === ParallelTransport() inverse_retraction_methods = [PolarInverseRetraction(), QRInverseRetraction()] diff --git a/test/manifolds/unitary_matrices.jl b/test/manifolds/unitary_matrices.jl index b0fc87cb06..ccd5a889f3 100644 --- a/test/manifolds/unitary_matrices.jl +++ b/test/manifolds/unitary_matrices.jl @@ -19,6 +19,7 @@ using Quaternions @test is_point(M, rand(MersenneTwister(), M); error=:error) @test abs(rand(MersenneTwister(), OrthogonalMatrices(1))[]) == 1 @test is_vector(M, p, rand(MersenneTwister(), M; vector_at=p)) + @test default_vector_transport_method(M) === ProjectionTransport() end @testset "Unitary Matrices" begin From 8d184415a4a94119dcd27bd92b6d6bc4706fe6b4 Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Sat, 23 Mar 2024 16:37:32 +0100 Subject: [PATCH 14/14] fix news --- NEWS.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index b806626109..faf84201c4 100644 --- a/NEWS.md +++ b/NEWS.md @@ -9,6 +9,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added +* using `DocumenterInterLinks` for links to other Julia packages documentation. * Implementation of `sectional_curvature`, `sectional_curvature_min` and `sectional_curvature_max` for several manifolds. * `sectional_curvature_matrix` function and a tutorial on coordinate-free curvature. @@ -18,7 +19,6 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Fixed -* using `DocumenterInterLinks` for links to other Julia packages documentation. * typographical errors in tutorials/working-in-charts.jl. * several typographical errors in the docs * unifies to use two backticks ``` `` ``` for math instead of ` $ ` further in the docs