🍩 (#533)

* Curved torus * formatting * general improvements * some work on tests and geodesics * partially working chart stitching * geodesic visualization * small fix * cleanup * faster solver * parallel transport * fix * cleanup * move test down * improve test coverage * polishing torus * small restructuring * log/distance using BVP in a chart * geodesic on torus in Makie * add some docs * Apply suggestions from code review Co-authored-by: Ronny Bergmann <[email protected]> * Addressing some review comments * addressing review * improve docs * bump version * minor improvements Co-authored-by: Ronny Bergmann <[email protected]>
JuliaManifolds · Oct 6, 2022 · e3af13e · e3af13e · mateuszbaran · Oct 6, 2022
1 parent c5b4990
commit e3af13e
Show file tree

Hide file tree

Showing 14 changed files with 977 additions and 6 deletions.
diff --git a/Project.toml b/Project.toml
@@ -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.8.30"
+version = "0.8.31"
 
 [deps]
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
@@ -32,7 +32,7 @@ Einsum = "0.4"
 Graphs = "1.4"
 HybridArrays = "0.4"
 Kronecker = "0.4, 0.5"
-ManifoldsBase = "0.13.17"
+ManifoldsBase = "0.13.19"
 MatrixEquations = "2.2"
 Plots = "1"
 Quaternions = "0.5"
@@ -46,7 +46,9 @@ StatsBase = "0.32, 0.33"
 julia = "1.6"
 
 [extras]
+BoundaryValueDiffEq = "764a87c0-6b3e-53db-9096-fe964310641d"
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
+DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
 DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78"
 FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
 Gtk = "4c0ca9eb-093a-5379-98c5-f87ac0bbbf44"
@@ -63,4 +65,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 VisualRegressionTests = "34922c18-7c2a-561c-bac1-01e79b2c4c92"
 
 [targets]
-test = ["Test", "Colors", "DoubleFloats", "FiniteDifferences", "Gtk", "ImageIO", "ImageMagick", "OrdinaryDiffEq", "NLsolve", "Plots", "PyPlot", "Quaternions", "QuartzImageIO", "RecipesBase"]
+test = ["Test", "BoundaryValueDiffEq", "Colors", "DiffEqCallbacks", "DoubleFloats", "FiniteDifferences", "Gtk", "ImageIO", "ImageMagick", "OrdinaryDiffEq", "NLsolve", "Plots", "PyPlot", "Quaternions", "QuartzImageIO", "RecipesBase"]
diff --git a/docs/Project.toml b/docs/Project.toml
@@ -1,23 +1,29 @@
 [deps]
+BoundaryValueDiffEq = "764a87c0-6b3e-53db-9096-fe964310641d"
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
+DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 HybridArrays = "1baab800-613f-4b0a-84e4-9cd3431bfbb9"
 Manifolds = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e"
 ManifoldsBase = "3362f125-f0bb-47a3-aa74-596ffd7ef2fb"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
+BoundaryValueDiffEq = "2"
+DiffEqCallbacks = "2"
 Documenter = "0.27"
 FiniteDifferences = "0.12"
 Graphs = "1.4"
 HybridArrays = "0.4"
 ManifoldsBase = "0.13"
+OrdinaryDiffEq = "6"
 Plots = "1"
 PyPlot = "2.9"
 StaticArrays = "1.0"
diff --git a/docs/make.jl b/docs/make.jl
@@ -1,6 +1,9 @@
 using Plots, RecipesBase, Manifolds, ManifoldsBase, Documenter, PyPlot
-# required for loading the manifold tests functios
+# required for loading methods that handle differential equation solving
+using OrdinaryDiffEq, BoundaryValueDiffEq, DiffEqCallbacks
+# required for loading the manifold tests functions
 using Test, FiniteDifferences
+
 ENV["GKSwstype"] = "100"
 
 generated_path = joinpath(@__DIR__, "src", "misc")

diff --git a/docs/src/features/atlases.md b/docs/src/features/atlases.md
@@ -49,3 +49,17 @@ Modules = [Manifolds,ManifoldsBase]
 Pages = ["cotangent_space.jl"]
 Order = [:type, :function]
 ```
+
+## Computations in charts
+
+```@autodocs
+Modules = [Manifolds,ManifoldsBase]
+Pages = ["ode_callback.jl"]
+Order = [:constant, :type, :function]
+```
+
+```@autodocs
+Modules = [Manifolds,ManifoldsBase]
+Pages = ["bvp.jl"]
+Order = [:constant, :type, :function]
+```
diff --git a/docs/src/manifolds/torus.md b/docs/src/manifolds/torus.md
@@ -24,3 +24,13 @@ Modules = [Manifolds]
 Pages = ["manifolds/Torus.jl"]
 Order = [:type, :function]
 ```
+
+## Embedded Torus
+
+Two-dimensional torus embedded in $ℝ^3$.
+
+```@autodocs
+Modules = [Manifolds]
+Pages = ["manifolds/EmbeddedTorus.jl"]
+Order = [:type, :function]
+```
diff --git a/src/Manifolds.jl b/src/Manifolds.jl
@@ -19,6 +19,8 @@ import ManifoldsBase:
     allocation_promotion_function,
     array_value,
     base_manifold,
+    change_basis,
+    change_basis!,
     check_point,
     check_size,
     check_vector,
@@ -325,6 +327,7 @@ include("manifolds/CenteredMatrices.jl")
 include("manifolds/CholeskySpace.jl")
 include("manifolds/Circle.jl")
 include("manifolds/Elliptope.jl")
+include("manifolds/EmbeddedTorus.jl")
 include("manifolds/FixedRankMatrices.jl")
 include("manifolds/GeneralizedGrassmann.jl")
 include("manifolds/GeneralizedStiefel.jl")
@@ -423,6 +426,11 @@ function Base.in(X, TpM::TangentSpaceAtPoint; kwargs...)
 end
 
 function __init__()
+    @require BoundaryValueDiffEq = "764a87c0-6b3e-53db-9096-fe964310641d" begin
+        using .BoundaryValueDiffEq
+        include("differentiation/bvp.jl")
+    end
+
     @require FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000" begin
         using .FiniteDifferences
         include("differentiation/finite_differences.jl")
@@ -431,6 +439,11 @@ function __init__()
     @require OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" begin
         using .OrdinaryDiffEq: ODEProblem, AutoVern9, Rodas5, solve
         include("differentiation/ode.jl")
+
+        @require DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def" begin
+            using .DiffEqCallbacks
+            include("differentiation/ode_callback.jl")
+        end
     end
 
     @require NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" begin
@@ -622,10 +635,13 @@ export ×,
     bundle_projection,
     canonical_project,
     canonical_project!,
+    change_basis,
+    change_basis!,
     change_metric,
     change_metric!,
     change_representer,
     change_representer!,
+    check_chart_switch,
     check_point,
     check_vector,
     christoffel_symbols_first,

diff --git a/src/atlases.jl b/src/atlases.jl
@@ -55,6 +55,26 @@ function RetractionAtlas(
 end
 RetractionAtlas() = RetractionAtlas(ExponentialRetraction(), LogarithmicInverseRetraction())
 
+"""
+    affine_connection(M::AbstractManifold, A::AbstractAtlas, i, a, Xc, Yc)
+
+Calculate affine connection on manifold `M` at point with parameters `a` in chart `i` of
+[`AbstractAtlas`](@ref) `A` of vectors with coefficients `Xc` and `Yc` in induced basis.
+"""
+function affine_connection(M::AbstractManifold, A, i, a, Xc, Yc)
+    Zc = allocate(Xc)
+    return affine_connection!(M, Zc, A, i, a, Xc, Yc)
+end
+
+"""
+    affine_connection!(M::AbstractManifold, Zc, A::AbstractAtlas, i, a, Xc, Yc)
+
+Calculate affine connection on manifold `M` at point with parameters `a` in chart `i` of an
+an [`AbstractAtlas`](@ref) `A` of vectors with coefficients `Zc` and `Yc` in induced basis and save the result
+in `Zc`.
+"""
+affine_connection!(M::AbstractManifold, Zc, A::AbstractAtlas, i, a, Xc, Yc)
+
 """
     get_default_atlas(::AbstractManifold)
 
@@ -97,6 +117,16 @@ end
     p,
 )
 
+"""
+    check_chart_switch(M::AbstractManifold, A::AbstractAtlas, i, a)
+
+Determine whether chart should be switched when an operation in chart `i` from an [`AbstractAtlas`](@ref) `A`
+reaches parameters `a` in that chart.
+
+By default `false` is returned.
+"""
+check_chart_switch(M::AbstractManifold, A::AbstractAtlas, i, a) = false
+
 function get_parameters!(M::AbstractManifold, a, A::RetractionAtlas, i, p)
     return get_coordinates!(M, a, i, inverse_retract(M, i, p, A.invretr), A.basis)
 end
@@ -158,6 +188,36 @@ get_chart_index(::AbstractManifold, ::AbstractAtlas, ::Any)
 
 get_chart_index(::AbstractManifold, ::RetractionAtlas, p) = p
 
+"""
+    get_chart_index(M::AbstractManifold, A::AbstractAtlas, i, a)
+
+Select a chart from an [`AbstractAtlas`](@ref) `A` for manifold `M` that is suitable for
+representing the neighborhood of point with parametrization `a` in chart `i`. This selection
+should be deterministic, although different charts may be selected for arbitrarily close but
+distinct points.
+
+# See also
+
+[`get_default_atlas`](@ref)
+"""
+get_chart_index(::AbstractManifold, ::AbstractAtlas, ::Any, ::Any)
+
+"""
+    inner(M::AbstractManifold, A::AbstractAtlas, i, a, Xc, Yc)
+
+Calculate inner product on manifold `M` at point with parameters `a` in chart `i` of an
+atlas `A` of vectors with coefficients `Xc` and `Yc` in induced basis.
+"""
+inner(M::AbstractManifold, A::AbstractAtlas, i, a, Xc, Yc)
+
+"""
+    norm(M::AbstractManifold, A::AbstractAtlas, i, a, Xc)
+
+Calculate norm on manifold `M` at point with parameters `a` in chart `i` of an
+[`AbstractAtlas`](@ref) `A` of vector with coefficients `Xc` in induced basis.
+"""
+norm(M::AbstractManifold, A::AbstractAtlas, i, a, Xc) = sqrt(inner(M, A, i, a, Xc, Xc))
+
 @doc raw"""
     transition_map(M::AbstractManifold, A_from::AbstractAtlas, i_from, A_to::AbstractAtlas, i_to, a)
     transition_map(M::AbstractManifold, A::AbstractAtlas, i_from, i_to, a)
@@ -213,6 +273,40 @@ function transition_map!(M::AbstractManifold, y, A::AbstractAtlas, i_from, i_to,
     return transition_map!(M, y, A, i_from, A, i_to, a)
 end
 
+"""
+    transition_map_diff(M::AbstractManifold, A::AbstractAtlas, i_from, a, c, i_to)
+
+Compute differential of transition map from chart `i_from` to chart `i_to` from an
+[`AbstractAtlas`](@ref) `A` on manifold `M` at point with parameters `a` on tangent vector
+with coordinates `c` in the induced basis.
+"""
+function transition_map_diff(M::AbstractManifold, A::AbstractAtlas, i_from, a, c, i_to)
+    old_B = induced_basis(M, A, i_from)
+    new_B = induced_basis(M, A, i_to)
+    p_final = get_point(M, A, i_from, a)
+    return change_basis(M, p_final, c, old_B, new_B)
+end
+
+"""
+    transition_map_diff!(M::AbstractManifold, c_out, A::AbstractAtlas, i_from, a, c, i_to)
+
+Compute [`transition_map_diff`](@ref) on given arguments and save the result in `c_out`.
+"""
+function transition_map_diff!(
+    M::AbstractManifold,
+    c_out,
+    A::AbstractAtlas,
+    i_from,
+    a,
+    c_in,
+    i_to,
+)
+    old_B = induced_basis(M, A, i_from)
+    new_B = induced_basis(M, A, i_to)
+    p_final = get_point(M, A, i_from, a)
+    return change_basis!(M, c_out, p_final, c_in, old_B, new_B)
+end
+
 """
     induced_basis(M::AbstractManifold, A::AbstractAtlas, i, p, VST::VectorSpaceType)
 
@@ -286,20 +380,31 @@ struct InducedBasis{𝔽,VST<:VectorSpaceType,TA<:AbstractAtlas,TI} <: AbstractB
 end
 
 """
-    induced_basis(::AbstractManifold, A::AbstractAtlas, i, VST::VectorSpaceType)
+    induced_basis(::AbstractManifold, A::AbstractAtlas, i, VST::VectorSpaceType = TangentSpace)
 
 Get the basis induced by chart with index `i` from an [`AbstractAtlas`](@ref) `A` of vector
 space of type `vs`. Returns an object of type [`InducedBasis`](@ref).
+
+# See also
+
+[`VectorSpaceType`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/bases.html#ManifoldsBase.VectorSpaceType), [`AbstractBasis`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/bases.html#ManifoldsBase.AbstractBasis)
 """
 function induced_basis(
     ::AbstractManifold{𝔽},
     A::AbstractAtlas,
     i,
-    VST::VectorSpaceType,
+    VST::VectorSpaceType=TangentSpace,
 ) where {𝔽}
     return InducedBasis{𝔽,typeof(VST),typeof(A),typeof(i)}(VST, A, i)
 end
 
+"""
+    inverse_chart_injectivity_radius(M::AbstractManifold, A::AbstractAtlas, i)
+
+Injectivity radius of `get_point` for chart `i` from an [`AbstractAtlas`](@ref) `A` of a manifold `M`.
+"""
+inverse_chart_injectivity_radius(M::AbstractManifold, A::AbstractAtlas, i)
+
 function dual_basis(
     M::AbstractManifold{𝔽},
     ::Any,