Merge branch 'kellertuer/steifel-canonical-metric' of github.com:Juli…

…aManifolds/Manifolds.jl into kellertuer/steifel-canonical-metric

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.4.17"
+version = "0.4.19"
 
 [deps]
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
@@ -48,6 +48,7 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Gtk = "4c0ca9eb-093a-5379-98c5-f87ac0bbbf44"
 ImageIO = "82e4d734-157c-48bb-816b-45c225c6df19"
 ImageMagick = "6218d12a-5da1-5696-b52f-db25d2ecc6d1"
+NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
@@ -59,4 +60,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 VisualRegressionTests = "34922c18-7c2a-561c-bac1-01e79b2c4c92"
 
 [targets]
-test = ["Test", "Colors", "DoubleFloats", "FiniteDiff", "ForwardDiff", "Gtk", "ImageIO", "ImageMagick", "OrdinaryDiffEq", "Plots", "PyPlot", "Quaternions", "QuartzImageIO", "RecipesBase", "ReverseDiff", "VisualRegressionTests"]
+test = ["Test", "Colors", "DoubleFloats", "FiniteDiff", "ForwardDiff", "Gtk", "ImageIO", "ImageMagick", "OrdinaryDiffEq", "NLsolve", "Plots", "PyPlot", "Quaternions", "QuartzImageIO", "RecipesBase", "ReverseDiff", "VisualRegressionTests"]

docs/src/manifolds/group.md

@@ -44,6 +44,22 @@ Pages = ["groups/circle_group.jl"]
 Order = [:type, :function]
 ```
 
+### General linear group
+
+```@autodocs
+Modules = [Manifolds]
+Pages = ["groups/general_linear.jl"]
+Order = [:type, :function]
+```
+
+### Special linear group
+
+```@autodocs
+Modules = [Manifolds]
+Pages = ["groups/special_linear.jl"]
+Order = [:type, :function]
+```
+
 ### Special orthogonal group
 
 ```@autodocs

docs/src/manifolds/stiefel.md

@@ -19,13 +19,6 @@ Order = [:function]
 ```
 
 ## The canonical metric
-
-```@autodocs
-Modules = [Manifolds]
-Pages = ["manifolds/StiefelCanonicalMetric.jl"]
-Order = [:type]
-```
-
 Any ``X∈T_p\mathcal M``, ``p∈\mathcal M``, can be written as
 
 ```math
@@ -36,9 +29,15 @@ A ∈ ℝ^{p×p} \text{ skew-symmetric},
 B ∈ ℝ^{n×p} \text{ arbitrary.}
 ```
 
-In the Euclidean case, the elements from ``A`` are counted twice (i.e. weighted with a factor of 2).
+In the [`EuclideanMetric`](@ref), the elements from ``A`` are counted twice (i.e. weighted with a factor of 2).
 The canonical metric avoids this.
 
+```@autodocs
+Modules = [Manifolds]
+Pages = ["manifolds/StiefelCanonicalMetric.jl"]
+Order = [:type]
+```
+
 ```@autodocs
 Modules = [Manifolds]
 Pages = ["manifolds/StiefelCanonicalMetric.jl"]

docs/src/misc/internals.md

@@ -7,10 +7,16 @@ This page documents the internal types and methods of `Manifolds.jl`'s that migh
 ```@docs
 Manifolds.eigen_safe
 Manifolds.find_pv
+Manifolds.isnormal
 Manifolds.log_safe
+Manifolds.log_safe!
+Manifolds.mul!_safe
 Manifolds.nzsign
+Manifolds.realify
+Manifolds.realify!
 Manifolds.size_to_tuple
 Manifolds.select_from_tuple
+Manifolds.unrealify!
 Manifolds.usinc
 Manifolds.usinc_from_cos
 Manifolds.vec2skew!

src/Manifolds.jl

@@ -132,6 +132,8 @@ include("manifolds/VectorBundle.jl")
 include("distributions.jl")
 include("projected_distribution.jl")
 
+include("approx_inverse_retraction.jl")
+
 # It's included early to ensure visibility of `Identity`
 include("groups/group.jl")
 
@@ -185,6 +187,8 @@ include("groups/array_manifold.jl")
 include("groups/product_group.jl")
 include("groups/semidirect_product_group.jl")
 
+include("groups/general_linear.jl")
+include("groups/special_linear.jl")
 include("groups/translation_group.jl")
 include("groups/special_orthogonal.jl")
 include("groups/circle_group.jl")
@@ -233,6 +237,11 @@ function __init__()
         include("ode.jl")
     end
 
+    @require NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" begin
+        using .NLsolve: NLsolve
+        include("nlsolve.jl")
+    end
+
     @require Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" begin
         using .Test: Test
         include("tests/tests_general.jl")
@@ -350,6 +359,7 @@ export AbstractRetractionMethod,
     ProductRetraction,
     PowerRetraction
 export AbstractInverseRetractionMethod,
+    ApproximateInverseRetraction,
     LogarithmicInverseRetraction,
     QRInverseRetraction,
     PolarInverseRetraction,
@@ -475,6 +485,7 @@ export AbstractGroupAction,
     ActionDirection,
     AdditionOperation,
     CircleGroup,
+    GeneralLinear,
     GroupManifold,
     GroupOperationAction,
     Identity,
@@ -489,6 +500,7 @@ export AbstractGroupAction,
     RotationAction,
     SemidirectProductGroup,
     SpecialEuclidean,
+    SpecialLinear,
     SpecialOrthogonal,
     TranslationGroup,
     TranslationAction

src/approx_inverse_retraction.jl

@@ -0,0 +1,104 @@
+"""
+    ApproximateInverseRetraction <: AbstractInverseRetractionMethod
+
+An abstract type for representing approximate inverse retraction methods.
+"""
+abstract type ApproximateInverseRetraction <: AbstractInverseRetractionMethod end
+
+"""
+    NLsolveInverseRetraction{T<:AbstractRetractionMethod,TV,TK} <:
+        ApproximateInverseRetraction
+
+An inverse retraction method for approximating the inverse of a retraction using `NLsolve`.
+
+# Constructor
+
+    NLsolveInverseRetraction(
+        method::AbstractRetractionMethod[, X0];
+        project_tangent=false,
+        project_point=false,
+        nlsolve_kwargs...,
+    )
+
+Constructs an approximate inverse retraction for the retraction `method` with initial guess
+`X0`, defaulting to the zero vector. If `project_tangent` is `true`, then the tangent
+vector is projected before the retraction using `project`. If `project_point` is `true`,
+then the resulting point is projected after the retraction. `nlsolve_kwargs` are keyword
+arguments passed to `NLsolve.nlsolve`.
+"""
+struct NLsolveInverseRetraction{TR<:AbstractRetractionMethod,TV,TK} <:
+       ApproximateInverseRetraction
+    retraction::TR
+    X0::TV
+    project_tangent::Bool
+    project_point::Bool
+    nlsolve_kwargs::TK
+    function NLsolveInverseRetraction(m, X0, project_point, project_tangent, nlsolve_kwargs)
+        isdefined(Manifolds, :NLsolve) ||
+            @warn "To use NLsolveInverseRetraction, NLsolve must be loaded using `using NLsolve`."
+        return new{typeof(m),typeof(X0),typeof(nlsolve_kwargs)}(
+            m,
+            X0,
+            project_point,
+            project_tangent,
+            nlsolve_kwargs,
+        )
+    end
+end
+function NLsolveInverseRetraction(
+    m,
+    X0=nothing;
+    project_tangent::Bool=false,
+    project_point::Bool=false,
+    nlsolve_kwargs...,
+)
+    return NLsolveInverseRetraction(m, X0, project_point, project_tangent, nlsolve_kwargs)
+end
+
+@decorator_transparent_signature inverse_retract(
+    M::AbstractDecoratorManifold,
+    p,
+    q,
+    m::NLsolveInverseRetraction,
+)
+
+@decorator_transparent_signature inverse_retract!(
+    M::AbstractDecoratorManifold,
+    X,
+    p,
+    q,
+    m::NLsolveInverseRetraction,
+)
+
+@doc raw"""
+    inverse_retract(M, p, q method::NLsolveInverseRetraction; kwargs...)
+
+Approximate the inverse of the retraction specified by `method.retraction` from `p` with
+respect to `q` on the [`Manifold`](@ref) `M` using NLsolve. This inverse retraction is
+not guaranteed to succeed and probably will not unless `q` is close to `p` and the initial
+guess `X0` is close.
+
+If the solver fails to converge, an [`OutOfInjectivityRadiusError`](@ref) is raised.
+See [`NLsolveInverseRetraction`](@ref) for configurable parameters.
+"""
+inverse_retract(::Manifold, p, q, ::NLsolveInverseRetraction; kwargs...)
+
+function inverse_retract!(M::Manifold, X, p, q, method::NLsolveInverseRetraction; kwargs...)
+    X0 = method.X0 === nothing ? zero_tangent_vector(M, p) : method.X0
+    res = _inverse_retract_nlsolve(
+        M,
+        p,
+        q,
+        method.retraction,
+        X0,
+        method.project_tangent,
+        method.project_point,
+        method.nlsolve_kwargs;
+        kwargs...,
+    )
+    if !res.f_converged
+        @debug res
+        throw(OutOfInjectivityRadiusError())
+    end
+    return copyto!(X, res.zero)
+end

src/groups/circle_group.jl

@@ -79,4 +79,4 @@ function group_log(G::CircleGroup, q)
     end
 end
 
-group_log!(G::CircleGroup, X, q) = (X .= group_log(G, q))
+group_log!(G::CircleGroup, X::AbstractVector, q::AbstractVector) = (X .= group_log(G, q))