🍩

Project.toml

@@ -47,6 +47,7 @@ julia = "1.6"
 
 [extras]
 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 +64,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", "Colors", "DiffEqCallbacks", "DoubleFloats", "FiniteDifferences", "Gtk", "ImageIO", "ImageMagick", "OrdinaryDiffEq", "NLsolve", "Plots", "PyPlot", "Quaternions", "QuartzImageIO", "RecipesBase"]

docs/src/manifolds/generalizedtorus.md

@@ -0,0 +1,13 @@
+# Generalized Torus
+
+Example: geodesic on a torus
+
+```@example
+
+```
+
+```@autodocs
+Modules = [Manifolds]
+Pages = ["manifolds/GeneralizedTorus.jl"]
+Order = [:type, :function]
+```

src/Manifolds.jl

@@ -328,6 +328,7 @@ include("manifolds/Elliptope.jl")
 include("manifolds/FixedRankMatrices.jl")
 include("manifolds/GeneralizedGrassmann.jl")
 include("manifolds/GeneralizedStiefel.jl")
+include("manifolds/GeneralizedTorus.jl")
 include("manifolds/Hyperbolic.jl")
 include("manifolds/MultinomialDoublyStochastic.jl")
 include("manifolds/MultinomialSymmetric.jl")
@@ -431,6 +432,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

src/atlases.jl

@@ -97,6 +97,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 atlas `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 +168,20 @@ 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)
+
 @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)
@@ -286,20 +310,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 atlas `A` of manifold `M`.
+"""
+inverse_chart_injectivity_radius(M::AbstractManifold, A::AbstractAtlas, i)
+
 function dual_basis(
     M::AbstractManifold{𝔽},
     ::Any,

src/differentiation/ode.jl

@@ -45,3 +45,46 @@ function exp!(
     copyto!(M, q, solve_exp_ode(M, p, X; kwargs...))
     return q
 end
+
+function chart_exp_problem(u, params, t)
+    M = params[1]
+    B = params[2]
+    p = u.x[1]
+    dx = u.x[2]
+    ddx = -affine_connection(M, p, dx, dx, B)
+    return ArrayPartition(dx, ddx)
+end
+
+function solve_chart_exp_ode(
+    M::AbstractManifold,
+    p,
+    X,
+    A::AbstractAtlas,
+    i;
+    final_time=1.0,
+    solver=AutoVern9(Rodas5()),
+    kwargs...,
+)
+    u0 = ArrayPartition(p, X)
+    B = induced_basis(M, A, i, TangentSpaceType())
+    params = (M, B)
+    prob = ODEProblem(chart_exp_problem, u0, (0.0, final_time), params)
+    sol = solve(prob, solver; kwargs...)
+    q = sol.u[end].x[1]
+    return sol
+    #println(sol)
+    return q
+end
+
+# also define exp! for metric manifold anew in this case
+function exp!(
+    ::TraitList{IsMetricManifold},
+    M::AbstractDecoratorManifold,
+    q,
+    p,
+    X;
+    kwargs...,
+)
+    copyto!(M, q, solve_exp_ode(M, p, X; kwargs...))
+    return q
+end

src/differentiation/ode_callback.jl

@@ -0,0 +1,85 @@
+
+"""
+    maybe_switch_chart(u, t, integrator)
+
+Terminate integration when integrator goes too closely to chart boundary.
+"""
+function maybe_switch_chart(u, t, integrator)
+    (M, B) = integrator.p
+    if check_chart_switch(M, B.A, B.i, u.x[1])
+        # switch charts
+        OrdinaryDiffEq.terminate!(integrator)
+    end
+    return u
+end
+
+"""
+    StitchedChartSolution{TM<:AbstractManifold,TA<:AbstractAtlas,TChart}
+
+Solution of an ODE on manifold `M` in charts of atlas `A`.
+"""
+struct StitchedChartSolution{TM<:AbstractManifold,TA<:AbstractAtlas,TChart}
+    M::TM
+    A::TA
+    sols::Vector{Tuple{OrdinaryDiffEq.SciMLBase.AbstractODESolution,TChart}}
+end
+
+function StitchedChartSolution(M::AbstractManifold, A::AbstractAtlas, TChart)
+    return StitchedChartSolution{typeof(M),typeof(A),TChart}(M, A, [])
+end
+
+function (scs::StitchedChartSolution)(t::Real)
+    for (sol, i) in scs.sols
+        if t <= sol.t[end]
+            B = induced_basis(scs.M, scs.A, i)
+            solt = sol(t)
+            p = get_point(scs.M, scs.A, i, solt.x[1])
+            X = get_vector(scs.M, p, solt.x[2], B)
+            return (p, X)
+        end
+    end
+    throw(DomainError("Time $t is outside of the solution."))
+end
+
+function (scs::StitchedChartSolution)(t::AbstractArray)
+    return map(scs, t)
+end
+
+function solve_chart_exp_ode(
+    M::AbstractManifold,
+    p,
+    X,
+    A::AbstractAtlas,
+    i0;
+    solver=AutoVern9(Rodas5()),
+    final_time=1.0,
+    kwargs...,
+)
+    u0 = ArrayPartition(p, X)
+    cur_i = i0
+    # callback stops solver when we get too close to chart boundary
+    cb = FunctionCallingCallback(maybe_switch_chart; func_start=false)
+    retcode = :Terminated
+    init_time = 0.0
+    sols = StitchedChartSolution(M, A, typeof(i0))
+    while retcode === :Terminated && init_time < final_time
+        B = induced_basis(M, A, cur_i)
+        params = (M, B)
+        prob =
+            ODEProblem(chart_exp_problem, u0, (init_time, final_time), params; callback=cb)
+        sol = solve(prob, solver; kwargs...)
+        retcode = sol.retcode
+        init_time = sol.t[end]
+        push!(sols.sols, (sol, cur_i))
+        # here we switch charts
+        new_i = get_chart_index(M, A, cur_i, sol.u[end].x[1])
+        new_p0 = transition_map(M, A, cur_i, new_i, sol.u[end].x[1])
+        new_B = induced_basis(M, A, new_i)
+        p_final = get_point(M, A, cur_i, sol.u[end].x[1])
+        new_X0 =
+            get_coordinates(M, p_final, get_vector(M, p_final, sol.u[end].x[2], B), new_B)
+        u0 = ArrayPartition(new_p0, new_X0)
+        cur_i = new_i
+    end
+    return sols
+end

src/manifolds/ConnectionManifold.jl

@@ -66,6 +66,32 @@ function active_traits(f, M::ConnectionManifold, args...)
     )
 end
 
+"""
+    affine_connection(M::AbstractManifold, p, X, Y)
+
+Calculate affine connection on manifold `M` at point `p` of vectors `X` and `Y`.
+"""
+function affine_connection(M::AbstractManifold, p, X, Y) end
+
+"""
+    affine_connection(M::AbstractManifold, p, Xc, Yc, B::AbstractBasis)
+
+Calculate affine connection on manifold `M` at point `p` of vectors with coefficients `Xc`
+and `Yc` in basis `B`.
+"""
+function affine_connection(M::AbstractManifold, p, Xc, Yc, B::AbstractBasis)
+    Zc = allocate(Xc)
+    return affine_connection!(M, Zc, p, Xc, Yc, B)
+end
+
+"""
+    affine_connection!(M::AbstractManifold, Zc, p, Xc, Yc, B::AbstractBasis)
+
+Calculate affine connection on manifold `M` at point `p` of vectors with coefficients `Xc`
+and `Yc` in basis `B` and save the result in `Zc`.
+"""
+function affine_connection!(M::AbstractManifold, Zc, p, Xc, Yc, B::AbstractBasis) end
+
 @doc raw"""
     christoffel_symbols_first(
         M::AbstractManifold,