move tensor product back to Manifolds.jl. move change_metric and repr…

…esenter from Manifolds.jl to here for the power manifold; reorganize test manifolds a bit

src/PowerManifold.jl

@@ -225,6 +225,47 @@ function allocation_promotion_function(M::AbstractPowerManifold, f, args::Tuple)
     return allocation_promotion_function(M.manifold, f, args)
 end
 
+"""
+    change_representer(M::AbstractPowerManifold, ::AbstractMetric, p, X)
+
+Since the metric on a power manifold decouples, the change of a representer can be done elementwise
+"""
+change_representer(::AbstractPowerManifold, ::AbstractMetric, ::Any, ::Any)
+
+function change_representer!(M::AbstractPowerManifold, Y, G::AbstractMetric, p, X)
+    rep_size = representation_size(M.manifold)
+    for i in get_iterator(M)
+        change_representer!(
+            M.manifold,
+            _write(M, rep_size, Y, i),
+            G,
+            _read(M, rep_size, p, i),
+            _read(M, rep_size, X, i),
+        )
+    end
+    return Y
+end
+
+"""
+    change_metric(M::AbstractPowerManifold, ::AbstractMetric, p, X)
+
+Since the metric on a power manifold decouples, the change of metric can be done elementwise.
+"""
+change_metric(M::AbstractPowerManifold, ::AbstractMetric, ::Any, ::Any)
+
+function change_metric!(M::AbstractPowerManifold, Y, G::AbstractMetric, p, X)
+    rep_size = representation_size(M.manifold)
+    for i in get_iterator(M)
+        change_metric!(
+            M.manifold,
+            _write(M, rep_size, Y, i),
+            G,
+            _read(M, rep_size, p, i),
+            _read(M, rep_size, X, i),
+        )
+    end
+    return Y
+end
 
 """
     check_point(M::AbstractPowerManifold, p; kwargs...)
@@ -1301,6 +1342,10 @@ function Base.show(
     return nothing
 end
 
+function vector_bundle_transport(::VectorSpaceType, ::PowerManifold)
+    return ParallelTransport()
+end
+
 function vector_transport_direction!(
     M::AbstractPowerManifold,
     Y,
@@ -1456,6 +1501,29 @@ function Base.view(
     return view(p[I...], rep_size_to_colons(rep_size)...)
 end
 
+@doc raw"""
+    Y = Weingarten(M::AbstractPowerManifold, p, X, V)
+    Weingarten!(M::AbstractPowerManifold, Y, p, X, V)
+
+Since the metric decouples, also the computation of the Weingarten map
+``\mathcal W_p`` can be computed elementwise on the single elements of the [`PowerManifold`](@ref) `M`.
+"""
+Weingarten(::AbstractPowerManifold, p, X, V)
+
+function Weingarten!(M::AbstractPowerManifold, Y, p, X, V)
+    rep_size = representation_size(M.manifold)
+    for i in get_iterator(M)
+        Weingarten!(
+            M.manifold,
+            _write(M, rep_size, Y, i),
+            _read(M, rep_size, p, i),
+            _read(M, rep_size, X, i),
+            _read(M, rep_size, V, i),
+        )
+    end
+    return Y
+end
+
 @inline function _write(M::AbstractPowerManifold, rep_size::Tuple, x::AbstractArray, i::Int)
     return _write(M, rep_size, x, (i,))
 end

src/VectorBundle.jl

@@ -344,9 +344,6 @@ function representation_size(B::TangentBundleFibers)
     return representation_size(B.manifold)
 end
 
-function Base.show(io::IO, tpt::TensorProductType)
-    return print(io, "TensorProductType(", join(tpt.spaces, ", "), ")")
-end
 function Base.show(io::IO, vb::VectorBundle)
     return print(io, "VectorBundle($(vb.type.fiber), $(vb.manifold))")
 end
@@ -368,14 +365,6 @@ Determine the vector tranport used for [`exp`](@ref exp(::FiberBundle, ::Any...)
 """
 fiber_bundle_transport(::VectorSpaceType, ::AbstractManifold) = ParallelTransport()
 
-function vector_space_dimension(M::AbstractManifold, V::TensorProductType)
-    dim = 1
-    for space in V.spaces
-        dim *= fiber_dimension(M, space)
-    end
-    return dim
-end
-
 function vector_space_dimension(B::TangentBundleFibers)
     return manifold_dimension(B.manifold)
 end

src/VectorFiber.jl

@@ -1,16 +1,4 @@
 
-"""
-    TensorProductType(spaces::VectorSpaceType...)
-
-Vector space type corresponding to the tensor product of given vector space
-types.
-"""
-struct TensorProductType{TS<:Tuple} <: VectorSpaceType
-    spaces::TS
-end
-
-TensorProductType(spaces::VectorSpaceType...) = TensorProductType{typeof(spaces)}(spaces)
-
 """
     VectorSpaceFiberType{TVS<:VectorSpaceType}
 

test/power.jl

@@ -5,6 +5,8 @@ using StaticArrays
 using LinearAlgebra
 using Random
 
+include("test_manifolds.jl")
+
 power_array_wrapper(::Type{NestedPowerRepresentation}, ::Int) = identity
 power_array_wrapper(::Type{NestedReplacingPowerRepresentation}, i::Int) = SVector{i}
 
@@ -293,4 +295,23 @@ struct TestArrayRepresentation <: AbstractPowerRepresentation end
                   rand(MersenneTwister(123), N; vector_at = p)
         end
     end
+
+    @testset "metric conversion" begin
+        M = TestSPD(3)
+        N = PowerManifold(M, NestedPowerRepresentation(), 2)
+        e = EuclideanMetric()
+        p = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1]
+        q = [2.0 0.0 0.0; 0.0 2.0 0.0; 0.0 0.0 1]
+        P = [p, q]
+        X = [log(M, p, q), log(M, q, p)]
+        Y = change_metric(N, e, P, X)
+        Yc = [change_metric(M, e, p, log(M, p, q)), change_metric(M, e, q, log(M, q, p))]
+        @test norm(N, P, Y .- Yc) ≈ 0
+        Z = change_representer(N, e, P, X)
+        Zc = [
+            change_representer(M, e, p, log(M, p, q)),
+            change_representer(M, e, q, log(M, q, p)),
+        ]
+        @test norm(N, P, Z .- Zc) ≈ 0
+    end
 end

test/test_manifolds.jl

@@ -0,0 +1,226 @@
+using ManifoldsBase: ℝ, ℂ, DefaultManifold, RealNumbers, EuclideanMetric
+using LinearAlgebra
+
+# minimal sphere implementation for testing more complicated manifolds
+struct TestSphere{N,𝔽} <: AbstractManifold{𝔽} end
+TestSphere(N::Int, 𝔽 = ℝ) = TestSphere{N,𝔽}()
+
+ManifoldsBase.representation_size(::TestSphere{N}) where {N} = (N + 1,)
+
+function ManifoldsBase.change_representer!(
+    M::TestSphere,
+    Y,
+    ::ManifoldsBase.EuclideanMetric,
+    p,
+    X,
+)
+    return copyto!(M, Y, p, X)
+end
+
+function ManifoldsBase.change_metric!(
+    M::TestSphere,
+    Y,
+    ::ManifoldsBase.EuclideanMetric,
+    p,
+    X,
+)
+    return copyto!(M, Y, p, X)
+end
+
+function ManifoldsBase.check_point(M::TestSphere, p; kwargs...)
+    if !isapprox(norm(p), 1.0; kwargs...)
+        return DomainError(
+            norm(p),
+            "The point $(p) does not lie on the $(M) since its norm is not 1.",
+        )
+    end
+    return nothing
+end
+
+function ManifoldsBase.check_vector(M::TestSphere, p, X; kwargs...)
+    if !isapprox(abs(real(dot(p, X))), 0.0; kwargs...)
+        return DomainError(
+            abs(dot(p, X)),
+            "The vector $(X) is not a tangent vector to $(p) on $(M), since it is not orthogonal in the embedding.",
+        )
+    end
+    return nothing
+end
+
+function ManifoldsBase.exp!(M::TestSphere, q, p, X)
+    return exp!(M, q, p, X, one(number_eltype(X)))
+end
+function ManifoldsBase.exp!(::TestSphere, q, p, X, t::Number)
+    θ = abs(t) * norm(X)
+    if θ == 0
+        copyto!(q, p)
+    else
+        X_scale = t * sin(θ) / θ
+        q .= p .* cos(θ) .+ X .* X_scale
+    end
+    return q
+end
+
+function ManifoldsBase.get_basis_diagonalizing(
+    M::TestSphere{n},
+    p,
+    B::DiagonalizingOrthonormalBasis{ℝ},
+) where {n}
+    A = zeros(n + 1, n + 1)
+    A[1, :] = transpose(p)
+    A[2, :] = transpose(B.frame_direction)
+    V = nullspace(A)
+    κ = ones(n)
+    if !iszero(B.frame_direction)
+        # if we have a nonzero direction for the geodesic, add it and it gets curvature zero from the tensor
+        V = hcat(B.frame_direction / norm(M, p, B.frame_direction), V)
+        κ[1] = 0 # no curvature along the geodesic direction, if x!=y
+    end
+    T = typeof(similar(B.frame_direction))
+    Ξ = [convert(T, V[:, i]) for i in 1:n]
+    return CachedBasis(B, κ, Ξ)
+end
+
+@inline function nzsign(z, absz = abs(z))
+    psignz = z / absz
+    return ifelse(iszero(absz), one(psignz), psignz)
+end
+
+
+function ManifoldsBase.get_coordinates_orthonormal!(M::TestSphere, Y, p, X, ::RealNumbers)
+    n = manifold_dimension(M)
+    p1 = p[1]
+    cosθ = abs(p1)
+    λ = nzsign(p1, cosθ)
+    pend, Xend = view(p, 2:(n + 1)), view(X, 2:(n + 1))
+    factor = λ * X[1] / (1 + cosθ)
+    Y .= Xend .- pend .* factor
+    return Y
+end
+
+function ManifoldsBase.get_vector_orthonormal!(M::TestSphere, Y, p, X, ::RealNumbers)
+    n = manifold_dimension(M)
+    p1 = p[1]
+    cosθ = abs(p1)
+    λ = nzsign(p1, cosθ)
+    pend = view(p, 2:(n + 1))
+    pX = dot(pend, X)
+    factor = pX / (1 + cosθ)
+    Y[1] = -λ * pX
+    Y[2:(n + 1)] .= X .- pend .* factor
+    return Y
+end
+
+ManifoldsBase.inner(::TestSphere, p, X, Y) = dot(X, Y)
+
+ManifoldsBase.injectivity_radius(::TestSphere) = π
+
+function ManifoldsBase.inverse_retract_project!(M::TestSphere, X, p, q)
+    X .= q .- p
+    project!(M, X, p, X)
+    return X
+end
+
+ManifoldsBase.is_flat(M::TestSphere) = manifold_dimension(M) == 1
+
+function ManifoldsBase.log!(::TestSphere, X, p, q)
+    cosθ = clamp(real(dot(p, q)), -1, 1)
+    X .= q .- p .* cosθ
+    nrmX = norm(X)
+    if nrmX > 0
+        X .*= acos(cosθ) / nrmX
+    end
+    return X
+end
+
+ManifoldsBase.manifold_dimension(::TestSphere{N}) where {N} = N
+
+function ManifoldsBase.parallel_transport_to!(::TestSphere, Y, p, X, q)
+    m = p .+ q
+    mnorm2 = real(dot(m, m))
+    factor = 2 * real(dot(X, q)) / mnorm2
+    Y .= X .- m .* factor
+    return Y
+end
+
+function ManifoldsBase.project!(::TestSphere, q, p)
+    q .= p ./ norm(p)
+    return q
+end
+
+function ManifoldsBase.project!(::TestSphere, Y, p, X)
+    Y .= X .- p .* real(dot(p, X))
+    return Y
+end
+
+function Random.rand!(M::TestSphere, pX; vector_at = nothing, σ = one(eltype(pX)))
+    return rand!(Random.default_rng(), M, pX; vector_at = vector_at, σ = σ)
+end
+function Random.rand!(
+    rng::AbstractRNG,
+    M::TestSphere,
+    pX;
+    vector_at = nothing,
+    σ = one(eltype(pX)),
+)
+    if vector_at === nothing
+        project!(M, pX, randn(rng, eltype(pX), representation_size(M)))
+    else
+        n = σ * randn(rng, eltype(pX), size(pX)) # Gaussian in embedding
+        project!(M, pX, vector_at, n) #project to TpM (keeps Gaussianness)
+    end
+    return pX
+end
+
+function ManifoldsBase.riemann_tensor!(M::TestSphere, Xresult, p, X, Y, Z)
+    innerZX = inner(M, p, Z, X)
+    innerZY = inner(M, p, Z, Y)
+    Xresult .= innerZY .* X .- innerZX .* Y
+    return Xresult
+end
+
+# from Absil, Mahony, Trumpf, 2013 https://sites.uclouvain.be/absil/2013-01/Weingarten_07PA_techrep.pdf
+function ManifoldsBase.Weingarten!(::TestSphere, Y, p, X, V)
+    return Y .= -X * p' * V
+end
+
+# minimal implementation of SPD (for its nontrivial metric conversion)
+
+struct TestSPD <: AbstractManifold{ℝ}
+    n::Int
+end
+
+function ManifoldsBase.change_representer!(::TestSPD, Y, ::EuclideanMetric, p, X)
+    Y .= p * X * p
+    return Y
+end
+
+function ManifoldsBase.change_metric!(::TestSPD, Y, ::EuclideanMetric, p, X)
+    Y .= p * X
+    return Y
+end
+
+function ManifoldsBase.inner(::TestSPD, p, X, Y)
+    F = cholesky(Symmetric(convert(AbstractMatrix, p)))
+    return dot((F \ Symmetric(X)), (Symmetric(Y) / F))
+end
+
+function spd_sqrt_and_sqrt_inv(p::AbstractMatrix)
+    e = eigen(Symmetric(p))
+    U = e.vectors
+    S = max.(e.values, floatmin(eltype(e.values)))
+    Ssqrt = Diagonal(sqrt.(S))
+    SsqrtInv = Diagonal(1 ./ sqrt.(S))
+    return (Symmetric(U * Ssqrt * transpose(U)), Symmetric(U * SsqrtInv * transpose(U)))
+end
+
+function ManifoldsBase.log!(::TestSPD, X, p, q)
+    (p_sqrt, p_sqrt_inv) = spd_sqrt_and_sqrt_inv(p)
+    T = Symmetric(p_sqrt_inv * convert(AbstractMatrix, q) * p_sqrt_inv)
+    e2 = eigen(T)
+    Se = Diagonal(log.(max.(e2.values, eps())))
+    pUe = p_sqrt * e2.vectors
+    return mul!(X, pUe, Se * transpose(pUe))
+end
+
+ManifoldsBase.representation_size(M::TestSPD) = (M.n, M.n)