rand for GL and Heisenberg group (#610)

* rand for GL and Heisenberg group

* Some docs

* format

* Update src/groups/general_linear.jl

Co-authored-by: Ronny Bergmann <git@ronnybergmann.net>

* fix compose on semidirect product group

* bump version

* Upload coverage from all Julia versions (less chance of random upload error leading to incomplete coverage)

---------

Co-authored-by: Ronny Bergmann <git@ronnybergmann.net>

.github/workflows/ci.yml

@@ -32,4 +32,4 @@ jobs:
       - uses: codecov/codecov-action@v3
         with:
           fail_ci_if_error: false
-        if: ${{ matrix.julia-version == '1.8' && matrix.os =='ubuntu-latest' }}
+        if: ${{ matrix.os =='ubuntu-latest' }}

Project.toml

@@ -1,7 +1,7 @@
 name = "Manifolds"
 uuid = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e"
 authors = ["Seth Axen <seth.axen@gmail.com>", "Mateusz Baran <mateuszbaran89@gmail.com>", "Ronny Bergmann <manopt@ronnybergmann.net>", "Antoine Levitt <antoine.levitt@gmail.com>"]
-version = "0.8.62"
+version = "0.8.63"
 
 [deps]
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"

src/groups/general_linear.jl

@@ -246,6 +246,26 @@ project(::GeneralLinear, p, X) = X
 project!(::GeneralLinear, q, p) = copyto!(q, p)
 project!(::GeneralLinear, Y, p, X) = copyto!(Y, X)
 
+@doc raw"""
+    Random.rand(G::GeneralLinear; vector_at=nothing, kwargs...)
+
+If `vector_at` is `nothing`, return a random point on the [`GeneralLinear`](@ref) group `G`
+by using `rand` in the embedding.
+
+If `vector_at` is not `nothing`, return a random tangent vector from the tangent space of
+the point `vector_at` on the [`GeneralLinear`](@ref) by using by using `rand` in the embedding.
+"""
+rand(G::GeneralLinear; kwargs...)
+
+function Random.rand!(G::GeneralLinear, pX; kwargs...)
+    rand!(get_embedding(G), pX; kwargs...)
+    return pX
+end
+function Random.rand!(rng::AbstractRNG, G::GeneralLinear, pX; kwargs...)
+    rand!(rng, get_embedding(G), pX; kwargs...)
+    return pX
+end
+
 Base.show(io::IO, ::GeneralLinear{n,𝔽}) where {n,𝔽} = print(io, "GeneralLinear($n, $𝔽)")
 
 translate_diff(::GeneralLinear, p, q, X, ::LeftAction) = X

src/groups/heisenberg.jl

@@ -353,6 +353,43 @@ function project!(M::HeisenbergGroup{n}, Y, p, X) where {n}
     return Y
 end
 
+@doc raw"""
+    Random.rand(M::HeisenbergGroup; vector_at = nothing, σ::Real=1.0)
+
+If `vector_at` is `nothing`, return a random point on the [`HeisenbergGroup`](@ref) `M`
+by sampling elements of the first row and the last column from the normal distribution with
+mean 0 and standard deviation `σ`.
+
+If `vector_at` is not `nothing`, return a random tangent vector from the tangent space of
+the point `vector_at` on the [`HeisenbergGroup`](@ref) by using a normal distribution with
+mean 0 and standard deviation `σ`.
+"""
+rand(M::HeisenbergGroup; vector_at=nothing, σ::Real=1.0)
+
+function Random.rand!(
+    rng::AbstractRNG,
+    ::HeisenbergGroup{n},
+    pX;
+    σ::Real=one(eltype(pX)),
+    vector_at=nothing,
+) where {n}
+    if vector_at === nothing
+        copyto!(pX, I)
+        va = view(pX, 1, 2:(n + 2))
+        randn!(rng, va)
+        va .*= σ
+        vb = view(pX, 2:(n + 1), n + 2)
+        randn!(rng, vb)
+        vb .*= σ
+    else
+        fill!(pX, 0)
+        randn!(rng, view(pX, 1, 2:(n + 2)))
+        randn!(rng, view(pX, 2:(n + 1), n + 2))
+        pX .*= σ
+    end
+    return pX
+end
+
 Base.show(io::IO, ::HeisenbergGroup{n}) where {n} = print(io, "HeisenbergGroup($n)")
 
 translate_diff(::HeisenbergGroup, p, q, X, ::LeftAction) = X

src/groups/semidirect_product_group.jl

@@ -123,13 +123,15 @@ function _compose!(G::SemidirectProductGroup, x, p, q)
     M = base_manifold(G)
     N, H = M.manifolds
     A = G.op.action
+    x_tmp = allocate(x)
     np, hp = submanifold_components(G, p)
     nq, hq = submanifold_components(G, q)
-    nx, hx = submanifold_components(G, x)
+    nx, hx = submanifold_components(G, x_tmp)
     compose!(H, hx, hp, hq)
     nxtmp = apply(A, hp, nq)
     compose!(N, nx, np, nxtmp)
     @inbounds _padpoint!(G, x)
+    copyto!(x, x_tmp)
     return x
 end
 

test/groups/general_linear.jl

@@ -131,6 +131,8 @@ using NLsolve
                 basis_types_to_from=basis_types,
                 exp_log_atol_multiplier=1e7,
                 retraction_atol_multiplier=1e7,
+                test_rand_point=true,
+                test_rand_tvector=true,
             )
         end
     end

test/groups/heisenberg.jl

@@ -52,5 +52,7 @@ include("group_utils.jl")
         test_project_point=true,
         test_project_tangent=true,
         test_musical_isomorphisms=false,
+        test_rand_point=true,
+        test_rand_tvector=true,
     )
 end

test/groups/semidirect_product_group.jl

@@ -43,6 +43,11 @@ include("group_utils.jl")
     @test compose(G, pts[1], e) == pts[1]
     @test compose(G, e, e) === e
 
+    # test in-place composition
+    o1 = copy(pts[1])
+    compose!(G, o1, o1, pts[2])
+    @test isapprox(G, o1, compose(G, pts[1], pts[2]))
+
     eA = identity_element(G)
     @test isapprox(G, eA, e)
     @test isapprox(G, e, eA)