Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Sectional curvature #711

Merged
merged 15 commits into from
Mar 23, 2024
4 changes: 2 additions & 2 deletions .github/workflows/documenter.yml
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
5 changes: 5 additions & 0 deletions NEWS.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
4 changes: 2 additions & 2 deletions Project.toml
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "Manifolds"
uuid = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e"
authors = ["Seth Axen <[email protected]>", "Mateusz Baran <[email protected]>", "Ronny Bergmann <[email protected]>", "Antoine Levitt <[email protected]>"]
version = "0.9.14"
version = "0.9.15"

[deps]
Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
Expand Down Expand Up @@ -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"
Expand Down
1 change: 1 addition & 0 deletions docs/make.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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" => [
Expand Down
12 changes: 3 additions & 9 deletions docs/src/features/utilities.md
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
1 change: 1 addition & 0 deletions docs/src/misc/notation.md
Original file line number Diff line number Diff line change
Expand Up @@ -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`` | | |
Expand Down
10 changes: 10 additions & 0 deletions docs/src/references.bib
Original file line number Diff line number Diff line change
Expand Up @@ -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},
Expand Down
7 changes: 7 additions & 0 deletions src/Manifolds.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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!,
Expand Down
27 changes: 27 additions & 0 deletions src/manifolds/Euclidean.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)")
Expand Down
39 changes: 39 additions & 0 deletions src/manifolds/Hyperbolic.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
42 changes: 42 additions & 0 deletions src/manifolds/Sphere.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down
24 changes: 24 additions & 0 deletions src/manifolds/SymmetricPositiveDefiniteAffineInvariant.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down
80 changes: 80 additions & 0 deletions src/utils.jl
Original file line number Diff line number Diff line change
Expand Up @@ -277,3 +277,83 @@ _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)
mateuszbaran marked this conversation as resolved.
Show resolved Hide resolved
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
22 changes: 11 additions & 11 deletions test/manifolds/embedded_torus.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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()
Expand Down
Loading
Loading