Add distance gradient (#15)

* add a gradient, fix some doc links.
* resolve #13
* resolve #14
* Rename inverse_retract_diff_argument_fd_approx to differential_inverse_retract_argument_fd_approx
* add Number and formatting.
* A note on the naming scheme.
* bump version to 0.3 due to the renaming.
* fix markdown.

Co-authored-by: Mateusz Baran <mateuszbaran89@gmail.com>

Project.toml

@@ -1,7 +1,7 @@
 name = "ManifoldDiff"
 uuid = "af67fdf4-a580-4b9f-bbec-742ef357defd"
 authors = ["Seth Axen <seth.axen@gmail.com>", "Mateusz Baran <mateuszbaran89@gmail.com>", "Ronny Bergmann <manopt@ronnybergmann.net>"]
-version = "0.2.1"
+version = "0.3"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

docs/Project.toml

@@ -8,16 +8,11 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 ManifoldDiff = "af67fdf4-a580-4b9f-bbec-742ef357defd"
 Manifolds = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e"
 ManifoldsBase = "3362f125-f0bb-47a3-aa74-596ffd7ef2fb"
-Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
-PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
 ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
-StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [compat]
 Documenter = "0.27"
-ManifoldDiff = "0.2"
+ManifoldDiff = "0.3"
 ManifoldsBase = "0.13"
-Manifolds = "0.8"
-Plots = "1"
-PyPlot = "2.9"
+Manifolds = "0.8"

docs/make.jl

@@ -1,7 +1,5 @@
-using Plots, Manifolds, ManifoldsBase, ManifoldDiff, Documenter, PyPlot
-ENV["GKSwstype"] = "100"
-# required for loading the manifold tests functions
-using Test, ForwardDiff, ReverseDiff, FiniteDifferences, Zygote
+using ForwardDiff, ReverseDiff, FiniteDifferences, Zygote
+using Manifolds, ManifoldsBase, ManifoldDiff, Documenter
 
 makedocs(
  # for development, we disable prettyurls

docs/src/index.md

@@ -1,3 +1,17 @@
 # ManifoldDiff
 
 The package __ManifoldDiff__ aims to provide automatic calculation of Riemannian gradients of functions defined on manifolds. It builds upon [`Manifolds.jl`](https://github.com/JuliaManifolds/Manifolds.jl).
+
+## Naming scheme
+
+Providing a derivative, differential or gradient for a given function, this package adds that information to the function name.
+For example
+
+* `grad_f` for a gradient ``\operatorname{grad} f``
+* `differential_f` for ``Df`` (also called pushforward)
+* `differential_f_variable` if `f` has multiple variables / parameters, since a usual writing in math is ``f_x`` in this case
+* `adjoint_differential_f` for pullbacks
+* `adjoint_differential_f_variable` if `f` has multiple variables / parameters
+* `f_derivative` for ``f'``
+
+the scheme is not completely fixed but tries to follow the mathematical notation.

docs/src/library.md

@@ -2,24 +2,45 @@
 
 Documentation for `ManifoldDiff.jl`'s methods and types for finite differences and automatic differentiation.
 
-## Differentials
+## Derivatives
 
 ```@autodocs
 Modules = [ManifoldDiff]
-Pages = ["differentials.jl"]
+Pages = ["derivatives.jl"]
 Order = [:type, :function, :constant]
+Private = true
 ```
 
+## Differentials and their adjoints
+
 ```@autodocs
 Modules = [ManifoldDiff]
 Pages = ["adjoint_differentials.jl"]
 Order = [:type, :function, :constant]
+Private = true
+```
+
+```@autodocs
+Modules = [ManifoldDiff]
+Pages = ["differentials.jl"]
+Order = [:type, :function, :constant]
+Private = true
 ```
 
 ```@autodocs
 Modules = [ManifoldDiff]
 Pages = ["diagonalizing_projectors.jl"]
 Order = [:type, :function, :constant]
+Private = true
+```
+
+## Gradients
+
+```@autodocs
+Modules = [ManifoldDiff]
+Pages = ["gradients.jl"]
+Order = [:type, :function, :constant]
+Private = true
 ```
 
 ## Jacobi fields
@@ -72,18 +93,12 @@ Pages = ["finite_differences.jl"]
 Order = [:type, :function, :constant]
 ```
 
-### ReverseDiff.jl
-
-```@autodocs
-Modules = [ManifoldDiff]
-Pages = ["reverse_diff.jl"]
-Order = [:type, :function, :constant]
-```
-
-### Zygote.jl
+## Internal functions
 
 ```@autodocs
 Modules = [ManifoldDiff]
-Pages = ["zygote.jl"]
+Pages = ["ManifoldDiff.jl"]
 Order = [:type, :function, :constant]
+Private = true
+Public=false
 ```

src/ManifoldDiff.jl

@@ -14,6 +14,7 @@ using ManifoldsBase:
  TangentSpaceType,
  PowerManifoldNested,
  PowerManifoldNestedReplacing,
+ allocate_result,
  get_iterator,
  _write,
  _read
@@ -36,7 +37,7 @@ struct NoneDiffBackend <: AbstractDiffBackend end
  _derivative(f, t[, backend::AbstractDiffBackend])
 
 Compute the derivative of a callable `f` at time `t` computed using the given `backend`,
-an object of type [`AbstractDiffBackend`](@ref). If the backend is not explicitly
+an object of type [`AbstractDiffBackend`](@ref ManifoldDiff.AbstractDiffBackend). If the backend is not explicitly
 specified, it is obtained using the function [`default_differential_backend`](@ref).
 
 This function calculates plain Euclidean derivatives, for Riemannian differentiation see
@@ -64,7 +65,7 @@ end
  _gradient(f, p[, backend::AbstractDiffBackend])
 
 Compute the gradient of a callable `f` at point `p` computed using the given `backend`,
-an object of type [`AbstractDiffBackend`](@ref). If the backend is not explicitly
+an object of type [`AbstractDiffBackend`](@ref ManifoldDiff.AbstractDiffBackend). If the backend is not explicitly
 specified, it is obtained using the function [`default_differential_backend`](@ref).
 
 This function calculates plain Euclidean gradients, for Riemannian gradient calculation see
@@ -164,13 +165,16 @@ end
 
 include("diagonalizing_projectors.jl")
 
-include("differentials.jl")
 include("adjoint_differentials.jl")
+include("derivatives.jl")
+include("differentials.jl")
+include("gradients.jl")
 include("Jacobi_fields.jl")
 
 include("riemannian_diff.jl")
 include("embedded_diff.jl")
 
+
 function __init__()
  @require FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" begin
  using .FiniteDiff

src/derivatives.jl

@@ -0,0 +1,68 @@
+@doc raw"""
+ Y = geodesic_derivative(M, p, X, t::Number; γt = geodesic(M, p, X, t))
+ geodesic_derivative!(M, Y, p, X, t::Number; γt = geodesic(M, p, X, t))
+
+Evaluate the derivative of the geodesic ``γ(t)`` with ``γ_{p,X}(0) = p`` and ``\dot γ_{p,X}(0) = X`` at ``t``.
+The formula reads
+
+```math
+\dot γ(t) = \mathcal P_{γ(t) \gets p} X
+```
+
+where $\mathcal P$ denotes the parallel transport.
+This computation can also be done in-place of ``Y``.
+
+# Optional Parameters
+
+* `γt` – (`geodesic(M, p, X, t)`) the point on the geodesic at ``t``.
+ This way if the point was computed earlier it can be resued here.
+"""
+function geodesic_derivative(M, p, X, t::Number; q = geodesic(M, p, X, t))
+ return parallel_transport_to(M, p, X, q)
+end
+function geodesic_derivative!(M, Y, p, X, t::Number; q = geodesic(M, p, X, t))
+ parallel_transport_to!(M, Y, p, X, q)
+ return Y
+end
+
+@doc raw"""
+ Y = shortest_geodesic_derivative(M, p, X, t::Number; γt = shortest_geodesic(M, p, q, t))
+ shortest_geodesic_derivative!(M, Y, p, X, t::Number; γt = shortest_geodesic(M, p, q, t))
+
+Evaluate the derivative of the shortest geodesic ``γ(t)`` with ``γ_{p,q}(0) = p`` and ``\dot γ_{p,q}(1) = q``
+at ``t``.
+The formula reads
+
+```math
+\dot γ(t) = \mathcal P_{γ(t) \gets p} \log_pq
+```
+
+where $\mathcal P$ denotes the parallel transport.
+This computation can also be done in-place of ``Y``.
+
+# Optional Parameters
+
+* `γt = geodesic(M, p, X, t)` the point on the geodesic at ``t``.
+ This way if the point was computed earlier it can be resued here.
+"""
+function shortest_geodesic_derivative(
+ M,
+ p,
+ q,
+ t::Number;
+ γt = shortest_geodesic(M, p, q, t),
+)
+ return parallel_transport_to(M, p, log(M, p, q), γt)
+end
+function shortest_geodesic_derivative!(
+ M,
+ Y,
+ p,
+ q,
+ t::Number;
+ γt = shortest_geodesic(M, p, q, t),
+)
+ log!(M, Y, p, q)
+ parallel_transport_to!(M, Y, p, Y, γt)
+ return Y
+end

src/diagonalizing_projectors.jl

@@ -18,7 +18,7 @@ Compute eigenvalues of the Jacobi operator $Y → R(Y,X)X$, where $R$ is the cur
 endomorphism, together with projectors onto eigenspaces of the operator.
 Projectors are objects of subtypes of [`AbstractProjector`](@ref).
 
-By default constructs projectors using the [`DiagonalizingOrthonormalBasis`](@ref).
+By default constructs projectors using the [`DiagonalizingOrthonormalBasis`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/bases.html#ManifoldsBase.DiagonalizingOrthonormalBasis).
 """
 function diagonalizing_projectors(M::AbstractManifold, p, X)
  B = get_basis(M, p, DiagonalizingOrthonormalBasis(X))
@@ -31,7 +31,7 @@ end
 
 """
  ProjectorOntoVector{TM<:AbstractManifold,TP,TX}
- 
+
 A structure that represents projector onto the subspace of the tangent space at `p`
 from manifold `M` spanned by tangent vector `X` of unit norm.
 
@@ -51,7 +51,7 @@ end
 
 """
  CoprojectorOntoVector{TM<:AbstractManifold,TP,TX}
- 
+
 
 A structure that represents projector onto the subspace of the tangent space at `p`
 from manifold `M` othogonal to vector `X` of unit norm.

src/differentials.jl

@@ -1,6 +1,5 @@
-
 @doc raw"""
- differential_shortest_geodesic_startpoint(M, p, q, t, X)
+ Y = differential_shortest_geodesic_startpoint(M, p, q, t, X)
  differential_shortest_geodesic_startpoint!(M, Y, p, q, t, X)
 
 Compute ``D_p γ(t;p,q)[η]`` (in place of `Y`).
@@ -17,7 +16,7 @@ end
 
 
 @doc raw"""
- differential_shortest_geodesic_endpoint(M, p, q, t, X)
+ Y = differential_shortest_geodesic_endpoint(M, p, q, t, X)
  differential_shortest_geodesic_endpoint!(M, Y, p, q, t, X)
 
 Compute ``D_qγ(t;p,q)[X]`` (in place of `Y`).
@@ -33,7 +32,7 @@ function differential_shortest_geodesic_endpoint!(M::AbstractManifold, Y, p, q,
 end
 
 @doc raw"""
- differential_exp_basepoint(M, p, X, Y)
+ Z = differential_exp_basepoint(M, p, X, Y)
  differential_exp_basepoint!(M, Z, p, X, Y)
 
 Compute ``D_p\exp_p X[Y]`` (in place of `Z`).
@@ -49,7 +48,7 @@ function differential_exp_basepoint!(M::AbstractManifold, Z, p, X, Y)
 end
 
 @doc raw"""
- differential_exp_argument(M, p, X, Y)
+ Z = differential_exp_argument(M, p, X, Y)
  differential_exp_argument!(M, Z, p, X, Y)
 
 computes ``D_X\exp_pX[Y]`` (in place of `Z`).
@@ -66,7 +65,7 @@ function differential_exp_argument!(M::AbstractManifold, Z, p, X, Y)
 end
 
 @doc raw"""
- differential_log_basepoint(M, p, q, X)
+ Y = differential_log_basepoint(M, p, q, X)
  differential_log_basepoint!(M, Y, p, q, X)
 
 computes ``D_p\log_pq[X]`` (in place of `Y`).
@@ -82,8 +81,8 @@ function differential_log_basepoint!(M::AbstractManifold, Y, p, q, X)
 end
 
 @doc raw"""
- differential_log_argument(M, p, q, X)
- differential_log_argument(M, Y, p, q, X)
+ Y = differential_log_argument(M, p, q, X)
+ differential_log_argument!(M, Y, p, q, X)
 
 computes ``D_q\log_pq[X]`` (in place of `Y`).
 
@@ -136,7 +135,7 @@ function differential_exp_argument_lie_approx!(M::AbstractManifold, Z, p, X, Y;
 end
 
 @doc raw"""
- inverse_retract_diff_argument_fd_approx(
+ differential_inverse_retract_argument_fd_approx(
  M::AbstractManifold,
  p,
  q,
@@ -159,7 +158,7 @@ retraction `invretr` and ``\operatorname{retr}`` is the retraction `retr`.
  > manifolds,” arXiv:1908.05875 [cs, math], Sep. 2019,
  > Available: http://arxiv.org/abs/1908.05875
 """
-function inverse_retract_diff_argument_fd_approx(
+function differential_inverse_retract_argument_fd_approx(
  M::AbstractManifold,
  p,
  q,
@@ -168,12 +167,12 @@ function inverse_retract_diff_argument_fd_approx(
  invretr::AbstractInverseRetractionMethod = default_inverse_retraction_method(M),
  h::Real = sqrt(eps(eltype(X))),
 )
- Y = allocate_result(M, inverse_retract_diff_argument_fd_approx, X, p, q)
- inverse_retract_diff_argument_fd_approx!(M, Y, p, q, X; retr, invretr, h)
+ Y = allocate_result(M, differential_inverse_retract_argument_fd_approx, X, p, q)
+ differential_inverse_retract_argument_fd_approx!(M, Y, p, q, X; retr, invretr, h)
  return Y
 end
 
-function inverse_retract_diff_argument_fd_approx!(
+function differential_inverse_retract_argument_fd_approx!(
  M::AbstractManifold,
  Y,
  p,