Fix random point / vector allocations (#653)

* fix allocation on tangent spaces and fixed rank for random point/vector generation.

* Fix positive numbers, bump dependencies.
* Fix SPDs
* fix Positive Numbers
* readd inplace rand for poisitive numbers and a corresponding test.
* bump version.

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.77"
+version = "0.8.78"
 
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
@@ -50,7 +50,7 @@ Graphs = "1.4"
 HybridArrays = "0.4"
 Kronecker = "0.4, 0.5"
 ManifoldDiff = "0.3.6"
-ManifoldsBase = "0.14.11"
+ManifoldsBase = "0.14.12"
 MatrixEquations = "2.2"
 OrdinaryDiffEq = "6.31"
 Plots = "1"

src/manifolds/FixedRankMatrices.jl

@@ -127,6 +127,17 @@ function allocate_result(
     # vals are p and X, so we can use their fields to set up those of the UMVTVector
     return UMVTVector(allocate(p.U, m, k), allocate(p.S, k, k), allocate(p.Vt, k, n))
 end
+function allocate_result(::FixedRankMatrices{m,n,k}, ::typeof(rand), p) where {m,n,k}
+    # vals are p and X, so we can use their fields to set up those of the UMVTVector
+    return UMVTVector(allocate(p.U, m, k), allocate(p.S, k, k), allocate(p.Vt, k, n))
+end
+function allocate_result(::FixedRankMatrices{m,n,k}, ::typeof(rand)) where {m,n,k}
+    return SVDMPoint(
+        Matrix{Float64}(undef, m, k),
+        Vector{Float64}(undef, k),
+        Matrix{Float64}(undef, k, n),
+    )
+end
 
 Base.copy(v::UMVTVector) = UMVTVector(copy(v.U), copy(v.M), copy(v.Vt))
 
@@ -448,32 +459,7 @@ and the singular values are sampled uniformly at random.
 If `vector_at` is not `nothing`, generate a random tangent vector in the tangent space of
 the point `vector_at` on the `FixedRankMatrices` manifold `M`.
 """
-function Random.rand(M::FixedRankMatrices; vector_at=nothing, kwargs...)
-    return rand(Random.default_rng(), M; vector_at=vector_at, kwargs...)
-end
-function Random.rand(
-    rng::AbstractRNG,
-    M::FixedRankMatrices{m,n,k};
-    vector_at=nothing,
-    kwargs...,
-) where {m,n,k}
-    if vector_at === nothing
-        p = SVDMPoint(
-            Matrix{Float64}(undef, m, k),
-            Vector{Float64}(undef, k),
-            Matrix{Float64}(undef, k, n),
-        )
-        return rand!(rng, M, p; kwargs...)
-    else
-        X = UMVTVector(
-            Matrix{Float64}(undef, m, k),
-            Matrix{Float64}(undef, k, k),
-            Matrix{Float64}(undef, k, n),
-        )
-        return rand!(rng, M, X; vector_at, kwargs...)
-    end
-end
-
+Random.rand(M::FixedRankMatrices; vector_at=nothing, kwargs...)
 function Random.rand!(
     rng::AbstractRNG,
     ::FixedRankMatrices{m,n,k},

src/manifolds/PositiveNumbers.jl

@@ -294,6 +294,19 @@ function parallel_transport_to!(::PositiveNumbers, Y, p, X, q)
     return (Y .= X .* q ./ p)
 end
 
+function Random.rand(M::PositiveNumbers; kwargs...)
+    return rand(Random.default_rng(), M; kwargs...)
+end
+
+function Random.rand(rng::AbstractRNG, ::PositiveNumbers; σ=1.0, vector_at=nothing)
+    if vector_at === nothing
+        pX = exp(randn(rng) * σ)
+    else
+        pX = vector_at * randn(rng) * σ
+    end
+    return pX
+end
+
 function Random.rand!(
     rng::AbstractRNG,
     ::PositiveNumbers,

src/manifolds/SymmetricPositiveDefinite.jl

@@ -403,6 +403,10 @@ Generate a random symmetric positive definite matrix on the
 """
 rand(M::SymmetricPositiveDefinite; σ::Real=1)
 
+function allocate_result(M::SymmetricPositiveDefinite, ::typeof(Random.rand), p::SPDPoint)
+    return zero_vector(M, p)
+end
+
 function Random.rand!(
     rng::AbstractRNG,
     M::SymmetricPositiveDefinite{N},

src/manifolds/VectorBundle.jl

@@ -292,6 +292,10 @@ See also [`VectorBundleInverseProductRetraction`](@ref).
 """
 struct VectorBundleProductRetraction <: AbstractRetractionMethod end
 
+function allocate_result(M::TangentSpaceAtPoint, ::typeof(rand))
+    return zero_vector(M.fiber.manifold, M.point)
+end
+
 base_manifold(B::VectorBundleFibers) = base_manifold(B.manifold)
 base_manifold(B::VectorSpaceAtPoint) = base_manifold(B.fiber)
 base_manifold(B::VectorBundle) = base_manifold(B.manifold)

test/manifolds/positive_numbers.jl

@@ -94,4 +94,12 @@ include("../utils.jl")
         @test volume_density(M, 0.5, 2.0) ≈ exp(4.0)
         @test volume_density(M5, [0.5, 1.0, 2.0], [1.0, -1.0, 2.0]) ≈ exp(2.0)
     end
+    @testset "Inplace random values" begin
+        p = fill(NaN)
+        X = fill(NaN)
+        rand!(M, p)
+        @test is_point(M, p)
+        rand!(M, X; vector_at=p)
+        @test is_vector(M, p, X)
+    end
 end