Add a numerical checks for retractions and their inverses (#187)

* Initial implementation, moving several methods over form Manopt * Unify loading and move. the last common type to the utils module. * Add documentation. * Update paper badge. * add vector check. * update News.md. * Apply suggestions from code review --------- Co-authored-by: Mateusz Baran <[email protected]>
JuliaManifolds · May 2, 2024 · 9d1af44 · 9d1af44
1 parent 830c5a1
commit 9d1af44
Show file tree

Hide file tree

Showing 23 changed files with 1,435 additions and 631 deletions.
diff --git a/NEWS.md b/NEWS.md
@@ -5,11 +5,18 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
-## [0.15.9] unreleased
+## [0.15.9] 2024-05-02
 
 ### Added
 
 * Tests now also use `Aqua.jl` to spot problems in the code such as ambiguities.
+* introduce a `check_inverse_retraction` function to numerically check whether an inverse retraction method is a (correct) inverse retraction.
+* introduce a `check_retraction` function to numerically check whether a retraction method is a (correct) retraction.
+* introduce a `check_vector_transport` function to numerically check whether a vector transport is a (correct) vector transport.
+
+### Changed
+
+* introduced a `ManifoldsBaseTestUtils` module to encapsulate common types and function definitions in different parts of the tests.
 
 ## [0.15.8] 13/03/2024
 

diff --git a/Project.toml b/Project.toml
@@ -1,44 +1,54 @@
 name = "ManifoldsBase"
 uuid = "3362f125-f0bb-47a3-aa74-596ffd7ef2fb"
 authors = ["Seth Axen <[email protected]>", "Mateusz Baran <[email protected]>", "Ronny Bergmann <[email protected]>", "Antoine Levitt <[email protected]>"]
-version = "0.15.8"
+version = "0.15.9"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 
 [weakdeps]
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [extensions]
+ManifoldsBasePlotsExt = "Plots"
 ManifoldsBaseRecursiveArrayToolsExt = "RecursiveArrayTools"
+ManifoldsBaseStatisticsExt = "Statistics"
 
 [compat]
 Aqua = "0.8"
 DoubleFloats = ">= 0.9.2"
 ForwardDiff = "0.10"
-julia = "1.6"
 LinearAlgebra = "1.6"
 Markdown = "1.6"
 OrdinaryDiffEq = "6"
+Plots = "1"
+Printf = "1.6"
 Random = "1.6"
 RecursiveArrayTools = "2, 3"
 Requires = "1"
 ReverseDiff = "1"
 StaticArrays = "1"
+Statistics = "1.6"
 Test = "1.6"
+julia = "1.6"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "Aqua", "DoubleFloats", "ForwardDiff", "OrdinaryDiffEq", "ReverseDiff", "StaticArrays", "RecursiveArrayTools"]
+test = ["Test", "Aqua", "DoubleFloats", "ForwardDiff", "OrdinaryDiffEq", "Plots", "ReverseDiff", "StaticArrays", "Statistics", "RecursiveArrayTools"]
diff --git a/Readme.md b/Readme.md
@@ -11,7 +11,7 @@
 [![codecov.io](http://codecov.io/github/JuliaManifolds/ManifoldsBase.jl/coverage.svg?branch=master)](https://codecov.io/gh/JuliaManifolds/ManifoldsBase.jl/)
 [![Aqua QA](https://raw.githubusercontent.com/JuliaTesting/Aqua.jl/master/badge.svg)](https://github.com/JuliaTesting/Aqua.jl)
 
-[![arXiv](https://img.shields.io/badge/arXiv%20CS.MS-2106.08777-blue.svg)](https://arxiv.org/abs/2106.08777)
+[![ACM TOMS](https://img.shields.io/badge/ACM%20TOMS-10.1145%2F3618296-blue.svg)](http://doi.org/10.1145/3618296)
 [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.5964340.svg)](https://doi.org/10.5281/zenodo.5964340)
 ## Installation
 

diff --git a/docs/make.jl b/docs/make.jl
@@ -82,6 +82,7 @@ makedocs(;
         "Meta-Manifolds" => "metamanifolds.md",
         "Decorating/Extending a Manifold" => "decorator.md",
         "Bases for tangent spaces" => "bases.md",
+        "Numerical Verification" => "numerical_verification.md",
         "References" => "references.md",
     ],
     plugins = [bib],

diff --git a/docs/src/numerical_verification.md b/docs/src/numerical_verification.md
@@ -0,0 +1,23 @@
+# Numerical Verification
+
+```@autodocs
+Modules = [ManifoldsBase]
+Pages = ["src/numerical_checks.jl"]
+Order = [:macro, :type, :function]
+Private = false
+Public = true
+```
+
+## Internal functions
+
+The following functions split the check into several parts, for
+example looking for the best fitting window and finding out the best slope,
+or plotting the slope.
+
+```@autodocs
+Modules = [ManifoldsBase]
+Pages = ["src/numerical_checks.jl"]
+Order = [:macro, :type, :function]
+Private = true
+Public = false
+```
diff --git a/docs/src/references.bib b/docs/src/references.bib
@@ -25,6 +25,18 @@ @article{AxenBaranBergmannRzecki:2023
     VOLUME    = {49},
     YEAR      = {2023}
 }
+@book{Boumal:2023,
+    TITLE     = {An Introduction to Optimization on Smooth Manifolds},
+    AUTHOR    = {Boumal, Nicolas},
+    YEAR      = {2023},
+    MONTH     = mar,
+    EDITION   = {First},
+    PUBLISHER = {Cambridge University Press},
+    DOI       = {10.1017/9781009166164},
+    URL       = {https://www.nicolasboumal.net/book/index.html},
+    ABSTRACT  = {Optimization on Riemannian manifolds-the result of smooth geometry and optimization merging into one elegant modern framework-spans many areas of science and engineering, including machine learning, computer vision, signal processing, dynamical systems and scientific computing. This text introduces the differential geometry and Riemannian geometry concepts that will help students and researchers in applied mathematics, computer science and engineering gain a firm mathematical grounding to use these tools confidently in their research. Its charts-last approach will prove more intuitive from an optimizer's viewpoint, and all definitions and theorems are motivated to build time-tested optimization algorithms. Starting from first principles, the text goes on to cover current research on topics including worst-case complexity and geodesic convexity. Readers will appreciate the tricks of the trade for conducting research and for numerical implementations sprinkled throughout the book.},
+    ISBN      = {978-1-00-916616-4}
+}
 @article{EhlersPiraniSchild:1972,
     DOI       = {10.1007/s10714-012-1353-4},
     YEAR      = {1972},

diff --git a/ext/ManifoldsBasePlotsExt.jl b/ext/ManifoldsBasePlotsExt.jl
@@ -0,0 +1,62 @@
+module ManifoldsBasePlotsExt
+
+if isdefined(Base, :get_extension)
+    using ManifoldsBase
+    using Plots
+    using Printf: @sprintf
+    import ManifoldsBase: plot_slope
+else
+    # imports need to be relative for Requires.jl-based workflows:
+    # https://github.com/JuliaArrays/ArrayInterface.jl/pull/387
+    using ..ManifoldsBase
+    using ..Plots
+    using ..ManifoldsBase: @sprintf # is from Printf, but loaded in ManifoldsBase, and since Printf loading works here only on full moon days between 12 and noon, this trick might do it?
+    import ..ManifoldsBase: plot_slope
+end
+
+function ManifoldsBase.plot_slope(
+    x,
+    y;
+    slope = 2,
+    line_base = 0,
+    a = 0,
+    b = 2.0,
+    i = 1,
+    j = length(x),
+)
+    fig = plot(
+        x,
+        y;
+        xaxis = :log,
+        yaxis = :log,
+        label = "\$E(t)\$",
+        linewidth = 3,
+        legend = :topleft,
+        color = :lightblue,
+    )
+    s_line = [exp10(line_base + t * slope) for t in log10.(x)]
+    plot!(
+        fig,
+        x,
+        s_line;
+        label = "slope s=$slope",
+        linestyle = :dash,
+        color = :black,
+        linewidth = 2,
+    )
+    if (i != 0) && (j != 0)
+        best_line = [exp10(a + t * b) for t in log10.(x[i:j])]
+        plot!(
+            fig,
+            x[i:j],
+            best_line;
+            label = "best slope $(@sprintf("%.4f", b))",
+            color = :blue,
+            linestyle = :dot,
+            linewidth = 2,
+        )
+    end
+    return fig
+end
+
+end
diff --git a/ext/ManifoldsBaseStatisticsExt.jl b/ext/ManifoldsBaseStatisticsExt.jl
@@ -0,0 +1,62 @@
+module ManifoldsBaseStatisticsExt
+
+if isdefined(Base, :get_extension)
+    using Statistics
+    using ManifoldsBase
+    import ManifoldsBase: find_best_slope_window
+else
+    # imports need to be relative for Requires.jl-based workflows:
+    # https://github.com/JuliaArrays/ArrayInterface.jl/pull/387
+    using ..Statistics
+    using ..ManifoldsBase
+    import ..ManifoldsBase: find_best_slope_window
+end
+
+function ManifoldsBase.find_best_slope_window(
+    X,
+    Y,
+    window = nothing;
+    slope::Real = 2.0,
+    slope_tol::Real = 0.1,
+)
+    n = length(X)
+    if window !== nothing && (any(window .> n))
+        error(
+            "One of the window sizes ($(window)) is larger than the length of the signal (n=$n).",
+        )
+    end
+    a_best = 0
+    b_best = -Inf
+    i_best = 0
+    j_best = 0
+    r_best = 0 # longest interval
+    for w in (window === nothing ? (2:n) : [window...])
+        for j in 1:(n - w + 1)
+            x = X[j:(j + w - 1)]
+            y = Y[j:(j + w - 1)]
+            # fit a line a + bx
+            c = cor(x, y)
+            b = std(y) / std(x) * c
+            a = mean(y) - b * mean(x)
+            # look for the largest interval where b is within slope tolerance
+            r = (maximum(x) - minimum(x))
+            if (r > r_best) && abs(b - slope) < slope_tol #longer interval found.
+                r_best = r
+                a_best = a
+                b_best = b
+                i_best = j
+                j_best = j + w - 1 #last index (see x and y from before)
+            end
+            # not best interval - maybe it is still the (first) best slope?
+            if r_best == 0 && abs(b - slope) < abs(b_best - slope)
+                # but do not update `r` since this indicates only a best r
+                a_best = a
+                b_best = b
+                i_best = j
+                j_best = j + w - 1 #last index (see x and y from before)
+            end
+        end
+    end
+    return (a_best, b_best, i_best, j_best)
+end
+end
diff --git a/src/ManifoldsBase.jl b/src/ManifoldsBase.jl
@@ -22,6 +22,7 @@ import Random: rand, rand!
 
 using LinearAlgebra
 using Markdown: @doc_str
+using Printf: @sprintf
 using Random
 using Requires
 
@@ -1056,6 +1057,7 @@ include("vector_spaces.jl")
 include("point_vector_fallbacks.jl")
 include("nested_trait.jl")
 include("decorator_trait.jl")
+include("numerical_checks.jl")
 
 include("VectorFiber.jl")
 include("TangentSpace.jl")
@@ -1067,15 +1069,53 @@ include("PowerManifold.jl")
 
 #
 #
-# Requires
+# Init
 # -----
 function __init__()
+    #
+    # Error Hints
+    #
+    @static if isdefined(Base.Experimental, :register_error_hint) # COV_EXCL_LINE
+        Base.Experimental.register_error_hint(MethodError) do io, exc, argtypes, kwargs
+            if exc.f === plot_slope
+                print(
+                    io,
+                    """
+
+                    `plot_slope` has to be implemented using your favourite plotting package.
+                    A default is available when Plots.jl is added to the current environment.
+                    To then get the plotting functionality activated, do
+                    """,
+                )
+                printstyled(io, "`using Plots`"; color = :cyan)
+            end
+            if exc.f === find_best_slope_window
+                print(
+                    io,
+                    """
+
+                    `find_best_slope_window` has to be implemented using some statistics package
+                    A default is available when Statistics.jl is added to the current environment.
+                    To then get the functionality activated, do
+                    """,
+                )
+                printstyled(io, "`using Statistics`"; color = :cyan)
+            end
+        end
+    end
+    # Extensions in the pre 1.9 fallback using Requires.jl
     @static if !isdefined(Base, :get_extension)
         @require RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" begin
             include(
                 "../ext/ManifoldsBaseRecursiveArrayToolsExt/ManifoldsBaseRecursiveArrayToolsExt.jl",
             )
         end
+        @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin
+            include("../ext/ManifoldsBasePlotsExt.jl")
+        end
+        @require Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" begin
+            include("../ext/ManifoldsBaseStatisticsExt.jl")
+        end
     end
 end
 #
@@ -1185,6 +1225,9 @@ export ×,
     change_metric!,
     change_representer,
     change_representer!,
+    check_inverse_retraction,
+    check_retraction,
+    check_vector_transport,
     copy,
     copyto!,
     default_approximation_method,