Introduce Kendall's shape space (#550)

* Introduce Kendall's shape space

* fixed shape space, added preshape space and some docs

* fix rotation action tests

* hand gesture example

* remove temporary `using Revise`

* improve tests and address an issue from review

* slightly expand tutorial

* more comments

* bump tolerance

* more docs

* bump tolerance to hopefully fix issue on Julia 1.6

* more reporting on why some tests fail

* increase tolerance

* addressing review

* more docs

* fix issues from CI

* improve coverage, bump version, fix incompatibility with newer OrdinaryDiffEq (termination codes changed)

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.38"
+version = "0.8.39"
 
 [deps]
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
@@ -34,6 +34,7 @@ HybridArrays = "0.4"
 Kronecker = "0.4, 0.5"
 ManifoldsBase = "0.13.21"
 MatrixEquations = "2.2"
+OrdinaryDiffEq = "6.31"
 Plots = "1"
 Quaternions = "0.5, 0.6"
 RecipesBase = "1.1"

docs/make.jl

@@ -99,6 +99,7 @@ makedocs(
                 "Projective space" => "manifolds/projectivespace.md",
                 "Orthogonal and Unitary Matrices" => "manifolds/generalunitary.md",
                 "Rotations" => "manifolds/rotations.md",
+                "Shape spaces" => "manifolds/shapespace.md",
                 "Skew-Hermitian matrices" => "manifolds/skewhermitian.md",
                 "Spectrahedron" => "manifolds/spectrahedron.md",
                 "Sphere" => "manifolds/sphere.md",

docs/src/manifolds/shapespace.md

@@ -0,0 +1,58 @@
+# Shape spaces
+
+Shape spaces are spaces of ``k`` points in ``\mathbb{R}^n`` up to simultaneous action of a group on all points.
+The most commonly encountered are Kendall's pre-shape and shape spaces.
+In the case of the Kendall's pre-shape spaces the action is translation and scaling.
+In the case of the Kendall's shape spaces the action is translation, scaling and rotation.
+
+```@example
+using Manifolds, Plots
+
+M = KendallsShapeSpace(2, 3)
+# two random point on the shape space
+p = [
+    0.4385117672460505 -0.6877826444042382 0.24927087715818771
+    -0.3830259932279294 0.35347460720654283 0.029551386021386548
+]
+q = [
+    -0.42693314765896473 -0.3268567431952937 0.7537898908542584
+    0.3054740561061169 -0.18962848284149897 -0.11584557326461796
+]
+# let's plot them as triples of points on a plane
+fig = scatter(p[1,:], p[2,:], label="p", aspect_ratio=:equal)
+scatter!(fig, q[1,:], q[2,:], label="q")
+
+# aligning q to p
+A = get_orbit_action(M)
+a = optimal_alignment(A, p, q)
+rot_q = apply(A, a, q)
+scatter!(fig, rot_q[1,:], rot_q[2,:], label="q aligned to p")
+```
+
+A more extensive usage example is available in the `hand_gestures.jl` tutorial.
+
+```@autodocs
+Modules = [Manifolds, ManifoldsBase]
+Pages = ["manifolds/KendallsPreShapeSpace.jl"]
+Order = [:type]
+```
+
+```@autodocs
+Modules = [Manifolds, ManifoldsBase]
+Pages = ["manifolds/KendallsShapeSpace.jl"]
+Order = [:type]
+```
+
+## Provided functions
+
+```@autodocs
+Modules = [Manifolds, ManifoldsBase]
+Pages = ["manifolds/KendallsPreShapeSpace.jl"]
+Order = [:function]
+```
+
+```@autodocs
+Modules = [Manifolds, ManifoldsBase]
+Pages = ["manifolds/KendallsShapeSpace.jl"]
+Order = [:function]
+```

src/Manifolds.jl

@@ -367,6 +367,10 @@ include("manifolds/Unitary.jl")
 include("manifolds/Rotations.jl")
 include("manifolds/Orthogonal.jl")
 
+# shape spaces require Sphere
+include("manifolds/KendallsPreShapeSpace.jl")
+include("manifolds/KendallsShapeSpace.jl")
+
 # Introduce the quotient, Grassmann, only after Stiefel
 include("manifolds/Grassmann.jl")
 
@@ -499,6 +503,8 @@ export Euclidean,
     Grassmann,
     HeisenbergGroup,
     Hyperbolic,
+    KendallsPreShapeSpace,
+    KendallsShapeSpace,
     Lorentz,
     MultinomialDoubleStochastic,
     MultinomialMatrices,

src/differentiation/ode_callback.jl

@@ -74,7 +74,11 @@ function (scs::StitchedChartSolution{:Exp})(t::Real)
             return (p, X)
         end
     end
-    throw(DomainError("Time $t is outside of the solution."))
+    throw(
+        DomainError(
+            "Time $t is outside of the solution (solution time range is [$(scs.sols[1][1].t[1]), $(scs.sols[end][1].t[end])]).",
+        ),
+    )
 end
 function (scs::StitchedChartSolution{:PT})(t::Real)
     if t < scs.sols[1][1].t[1]
@@ -90,7 +94,11 @@ function (scs::StitchedChartSolution{:PT})(t::Real)
             return (p, X, Y)
         end
     end
-    throw(DomainError("Time $t is outside of the solution."))
+    throw(
+        DomainError(
+            "Time $t is outside of the solution (solution time range is [$(scs.sols[1][1].t[1]), $(scs.sols[end][1].t[end])]).",
+        ),
+    )
 end
 
 function (scs::StitchedChartSolution)(t::AbstractArray)
@@ -139,10 +147,10 @@ function solve_chart_exp_ode(
         IntegratorTerminatorNearChartBoundary(check_chart_switch_kwargs);
         func_start=false,
     )
-    retcode = :Terminated
+    retcode = SciMLBase.ReturnCode.Terminated
     init_time = zero(final_time)
     sols = StitchedChartSolution(M, A, :Exp, typeof(i0))
-    while retcode === :Terminated && init_time < final_time
+    while retcode === SciMLBase.ReturnCode.Terminated && init_time < final_time
         params = (M, A, cur_i)
         prob =
             ODEProblem(chart_exp_problem, u0, (init_time, final_time), params; callback=cb)
@@ -216,10 +224,10 @@ function solve_chart_parallel_transport_ode(
         IntegratorTerminatorNearChartBoundary(check_chart_switch_kwargs);
         func_start=false,
     )
-    retcode = :Terminated
+    retcode = SciMLBase.ReturnCode.Terminated
     init_time = zero(final_time)
     sols = StitchedChartSolution(M, A, :PT, typeof(i0))
-    while retcode === :Terminated && init_time < final_time
+    while retcode === SciMLBase.ReturnCode.Terminated && init_time < final_time
         params = (M, A, cur_i)
         prob =
             ODEProblem(chart_pt_problem, u0, (init_time, final_time), params; callback=cb)

src/groups/rotation_action.jl

@@ -214,3 +214,100 @@ group_manifold(A::RowwiseMultiplicationAction) = A.manifold
 function inverse_apply(::LeftRowwiseMultiplicationAction, a, p)
     return (a \ p')'
 end
+
+###
+
+@doc raw"""
+    ColumnwiseMultiplicationAction{
+        TM<:AbstractManifold,
+        TO<:GeneralUnitaryMultiplicationGroup,
+        TAD<:ActionDirection,
+    } <: AbstractGroupAction{TAD}
+
+Action of the (special) unitary or orthogonal group [`GeneralUnitaryMultiplicationGroup`](@ref)
+of type `On` columns of points on a matrix manifold `M`.
+
+# Constructor
+
+    ColumnwiseMultiplicationAction(
+        M::AbstractManifold,
+        On::GeneralUnitaryMultiplicationGroup,
+        AD::ActionDirection = LeftAction(),
+    )
+"""
+struct ColumnwiseMultiplicationAction{
+    TM<:AbstractManifold,
+    TO<:GeneralUnitaryMultiplicationGroup,
+    TAD<:ActionDirection,
+} <: AbstractGroupAction{TAD}
+    manifold::TM
+    On::TO
+end
+
+function ColumnwiseMultiplicationAction(
+    M::AbstractManifold,
+    On::GeneralUnitaryMultiplicationGroup,
+    ::TAD=LeftAction(),
+) where {TAD<:ActionDirection}
+    return ColumnwiseMultiplicationAction{typeof(M),typeof(On),TAD}(M, On)
+end
+
+const LeftColumnwiseMultiplicationAction{
+    TM<:AbstractManifold,
+    TO<:GeneralUnitaryMultiplicationGroup,
+} = ColumnwiseMultiplicationAction{TM,TO,LeftAction}
+
+function apply(::LeftColumnwiseMultiplicationAction, a, p)
+    return a * p
+end
+function apply(::LeftColumnwiseMultiplicationAction, ::Identity{MultiplicationOperation}, p)
+    return p
+end
+
+function apply!(::LeftColumnwiseMultiplicationAction, q, a, p)
+    return map((qrow, prow) -> mul!(qrow, a, prow), eachcol(q), eachcol(p))
+end
+
+base_group(A::LeftColumnwiseMultiplicationAction) = A.On
+
+group_manifold(A::LeftColumnwiseMultiplicationAction) = A.manifold
+
+function inverse_apply(::LeftColumnwiseMultiplicationAction, a, p)
+    return a \ p
+end
+
+@doc raw"""
+    optimal_alignment(A::LeftColumnwiseMultiplicationAction, p, q)
+
+Compute optimal alignment for the left [`ColumnwiseMultiplicationAction`](@ref), i.e. the
+group element ``O^{*}`` that, when it acts on `p`, returns the point closest to `q`. Details
+of computation are described in Section 2.2.1 of [^Srivastava2016].
+
+The formula reads
+```math
+O^{*} = \begin{cases}
+UV^T & \text{if } \operatorname{det}(p q^{\mathrm{T}})\\
+U K V^{\mathrm{T}} & \text{otherwise}
+\end{cases}
+```
+where ``U \Sigma V^{\mathrm{T}}`` is the SVD decomposition of ``p q^{\mathrm{T}}`` and ``K``
+is the unit diagonal matrix with the last element on the diagonal replaced with -1.
+
+# References
+
+[^Srivastava2016]:
+    > A. Srivastava and E. P. Klassen, Functional and Shape Data Analysis. Springer New York, 2016.
+    > ISBN: 978-1-4939-4018-9.
+    > doi: [10.1007/978-1-4939-4020-2](https://doi.org/10.1007/978-1-4939-4020-2).
+"""
+function optimal_alignment(A::LeftColumnwiseMultiplicationAction, p, q)
+    is_point(A.manifold, p, true)
+    is_point(A.manifold, q, true)
+
+    Xmul = p * transpose(q)
+    F = svd(Xmul)
+    L = size(Xmul)[2]
+    UVt = F.U * F.Vt
+    Ostar = det(UVt) ≥ 0 ? UVt : F.U * Diagonal([i < L ? 1 : -1 for i in 1:L]) * F.Vt
+    return convert(typeof(Xmul), Ostar)
+end