diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 364688c3d4..9344289ffd 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -32,5 +32,8 @@ jobs: - uses: julia-actions/julia-processcoverage@v1 - uses: codecov/codecov-action@v4 with: + token: ${{ secrets.CODECOV_TOKEN }} + file: ./lcov.info + name: codecov-umbrella fail_ci_if_error: false if: ${{ matrix.os =='ubuntu-latest' }} diff --git a/NEWS.md b/NEWS.md index 14a4a549c6..23af11c3d3 100644 --- a/NEWS.md +++ b/NEWS.md @@ -5,7 +5,13 @@ 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.9.16] – unreleased +## [0.9.17] – unreleased + +### Added + +* `Hyperrectangle` manifold with boundary. + +## [0.9.16] – 2024-04-01 ### Changed diff --git a/Project.toml b/Project.toml index 814d794b7d..33b01c3a45 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Manifolds" uuid = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e" authors = ["Seth Axen ", "Mateusz Baran ", "Ronny Bergmann ", "Antoine Levitt "] -version = "0.9.16" +version = "0.9.17" [deps] Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" diff --git a/docs/make.jl b/docs/make.jl index 838774f93c..22713d86ff 100755 --- a/docs/make.jl +++ b/docs/make.jl @@ -121,6 +121,7 @@ makedocs(; "Grassmann" => "manifolds/grassmann.md", "Hamiltonian" => "manifolds/hamiltonian.md", "Hyperbolic space" => "manifolds/hyperbolic.md", + "Hyperrectangle" => "manifolds/hyperrectangle.md", "Lorentzian manifold" => "manifolds/lorentz.md", "Multinomial doubly stochastic matrices" => "manifolds/multinomialdoublystochastic.md", "Multinomial matrices" => "manifolds/multinomial.md", diff --git a/docs/src/manifolds/grassmann.md b/docs/src/manifolds/grassmann.md index c3dc8bfc28..2c4d7900a3 100644 --- a/docs/src/manifolds/grassmann.md +++ b/docs/src/manifolds/grassmann.md @@ -22,10 +22,9 @@ Pages = ["manifolds/GrassmannProjector.jl"] Order = [:type,:function] ``` - ## Literature ```@bibliography Pages = ["grassmann.md"] Canonical=false -``` \ No newline at end of file +``` diff --git a/docs/src/manifolds/hyperrectangle.md b/docs/src/manifolds/hyperrectangle.md new file mode 100644 index 0000000000..cad7119ad6 --- /dev/null +++ b/docs/src/manifolds/hyperrectangle.md @@ -0,0 +1,20 @@ +# [Hyperrectangle](@id HyperrectangleSection) + +Hyperrectangle is a manifold with corners [Joyce:2010](@cite), and also a subset of the real [Euclidean](@ref Main.Manifolds.Euclidean) manifold. +It is useful for box-constrained optimization, for example it is implicitly used in the classic L-BFGS-B algorithm. + +!!! note + This is a manifold with corners. Some parts of its interface specific to this property are experimental and may change without a breaking release. + +```@autodocs +Modules = [Manifolds] +Pages = ["manifolds/Hyperrectangle.jl"] +Order = [:type,:function] +``` + +## Literature + +```@bibliography +Pages = ["hyperrectangle.md"] +Canonical=false +``` diff --git a/docs/src/references.bib b/docs/src/references.bib index e377fd0737..8596c243c8 100644 --- a/docs/src/references.bib +++ b/docs/src/references.bib @@ -453,6 +453,15 @@ @article{JourneeBachAbsilSepulchre:2010 TITLE = {Low-Rank Optimization on the Cone of Positive Semidefinite Matrices}, JOURNAL = {SIAM Journal on Optimization} } + +@misc{Joyce:2010, + TITLE = {On manifolds with corners}, + DOI = {10.48550/arXiv.0910.3518}, + AUTHOR = {Joyce, Dominic}, + YEAR = {2010}, + EPRINT = {0910.3518}, + EPRINTTYPE = {arXiv} +} # # K # ---------------------------------------------------------------------------------------- diff --git a/src/Manifolds.jl b/src/Manifolds.jl index 317c663cef..f9af303499 100644 --- a/src/Manifolds.jl +++ b/src/Manifolds.jl @@ -428,6 +428,7 @@ include("manifolds/FlagStiefel.jl") include("manifolds/GeneralizedGrassmann.jl") include("manifolds/GeneralizedStiefel.jl") include("manifolds/Hyperbolic.jl") +include("manifolds/Hyperrectangle.jl") include("manifolds/MultinomialDoublyStochastic.jl") include("manifolds/MultinomialSymmetric.jl") include("manifolds/MultinomialSymmetricPositiveDefinite.jl") @@ -645,6 +646,7 @@ export Euclidean, HamiltonianMatrices, HeisenbergGroup, Hyperbolic, + Hyperrectangle, KendallsPreShapeSpace, KendallsShapeSpace, Lorentz, diff --git a/src/manifolds/Euclidean.jl b/src/manifolds/Euclidean.jl index 0098faa3a5..e2b471fad2 100644 --- a/src/manifolds/Euclidean.jl +++ b/src/manifolds/Euclidean.jl @@ -407,7 +407,7 @@ where ``I`` is the set of vectors ``k ∈ ℕ^i``, such that for all For the special case of ``i ≤ 2``, i.e. matrices and vectors, this simplifies to ````math -g_p(X,Y) = X^{\mathrm{H}}Y, +g_p(X,Y) = \operatorname{tr}(X^{\mathrm{H}}Y), ```` where ``⋅^{\mathrm{H}}`` denotes the Hermitian, i.e. complex conjugate transposed. @@ -745,7 +745,6 @@ end Compute the Riemann tensor ``R(X,Y)Z`` at point `p` on [`Euclidean`](@ref) manifold `M`. Its value is always the zero tangent vector. -```` """ riemann_tensor(M::Euclidean, p, X, Y, Z) @@ -894,13 +893,13 @@ Weingarten(::Euclidean, p, X, V) Weingarten!(::Euclidean, Y, p, X, V) = fill!(Y, 0) """ - zero_vector(M::Euclidean, x) + zero_vector(M::Euclidean, p) -Return the zero vector in the tangent space of `x` on the [`Euclidean`](@ref) -`M`, which here is just a zero filled array the same size as `x`. +Return the zero vector in the tangent space of `p` on the [`Euclidean`](@ref) +`M`, which here is just a zero filled array the same size as `p`. """ zero_vector(::Euclidean, ::Any...) zero_vector(::Euclidean{TypeParameter{Tuple{}}}, p::Number) = zero(p) zero_vector(::Euclidean{Tuple{}}, p::Number) = zero(p) -zero_vector!(::Euclidean, v, ::Any) = fill!(v, 0) +zero_vector!(::Euclidean, X, ::Any) = fill!(X, 0) diff --git a/src/manifolds/Hyperrectangle.jl b/src/manifolds/Hyperrectangle.jl new file mode 100644 index 0000000000..11c0180fae --- /dev/null +++ b/src/manifolds/Hyperrectangle.jl @@ -0,0 +1,530 @@ +@doc raw""" + Hyperrectangle{T} <: AbstractManifold{ℝ} + +Hyperrectangle, also known as orthotope or box. This is a manifold with corners [Joyce:2010](@cite) +with the standard Euclidean metric. + +# Constructor + + Hyperrectangle(lb::AbstractArray, ub::AbstractArray) + +Generate the hyperrectangle of arrays such that each element of the array is between lower +and upper bound with the same index. +""" +struct Hyperrectangle{T<:AbstractArray} <: AbstractDecoratorManifold{ℝ} + lb::T + ub::T + function Hyperrectangle(lb::T, ub::T) where {T<:AbstractArray} + for i in eachindex(lb, ub) + if lb[i] > ub[i] + throw( + ArgumentError( + "Lower bound at index $i ($(lb[i])) is greater than the upper bound ($(ub[i]))", + ), + ) + end + end + return new{T}(lb, ub) + end +end + +function active_traits(f, ::Hyperrectangle, args...) + return merge_traits( + IsDefaultMetric(EuclideanMetric()), + IsDefaultConnection(LeviCivitaConnection()), + ) +end + +function check_point(M::Hyperrectangle, p) + if !(eltype(p) <: Real) + return DomainError( + eltype(p), + "The matrix $(p) is not a real-valued matrix, so it does not lie on $(M).", + ) + end + for i in eachindex(M.lb, M.ub, p) + if p[i] < M.lb[i] + return DomainError( + p[i], + "At index $i the point has coordinate $(p[i]), below the lower bound of $(M.lb[i])", + ) + end + if p[i] > M.ub[i] + return DomainError( + p[i], + "At index $i the point has coordinate $(p[i]), above the upper bound of $(M.ub[i])", + ) + end + end + return nothing +end + +function check_vector(M::Hyperrectangle, p, X; kwargs...) + if !(eltype(X) <: Real) + return DomainError( + eltype(X), + "The matrix $(X) is not a real-valued matrix, so it can not be a tangent vector to $(p) on $(M).", + ) + end + return nothing +end + +default_approximation_method(::Hyperrectangle, ::typeof(mean)) = EfficientEstimator() + +""" + default_retraction_method(M::Hyperrectangle) + +Return [`ProjectionRetraction`](@extref `ManifoldsBase.ProjectionRetraction`) as the default +retraction for the [`Hyperrectangle`](@ref) manifold. +""" +default_retraction_method(::Hyperrectangle) = ProjectionRetraction() + +""" + distance(M::Hyperrectangle, p, q) + +Compute the euclidean distance between two points on the [`Hyperrectangle`](@ref) +manifold `M`, i.e. for vectors it's just the norm of the difference, for matrices +and higher order arrays, the matrix and tensor Frobenius norm, respectively. +""" +Base.@propagate_inbounds function distance(M::Hyperrectangle, p, q) + # Inspired by euclidean distance calculation in Distances.jl + # Much faster for large p, q than a naive implementation + @boundscheck if axes(p) != axes(q) + throw(DimensionMismatch("At last one of $p and $q does not belong to $M")) + end + s = zero(eltype(p)) + @inbounds begin + @simd for I in eachindex(p, q) + p_i = p[I] + q_i = q[I] + s += abs2(p_i - q_i) + end # COV_EXCL_LINE + end + return sqrt(s) +end +distance(::Hyperrectangle, p::Number, q::Number) = abs(p - q) + +""" + embed(M::Hyperrectangle, p) + +Embed the point `p` in `M`. Equivalent to an identity map. +""" +embed(::Hyperrectangle, p) = p + +""" + embed(M::Hyperrectangle, p, X) + +Embed the tangent vector `X` at point `p` in `M`. Equivalent to an identity map. +""" +embed(::Hyperrectangle, p, X) = X + +@doc raw""" + exp(M::Hyperrectangle, p, X) + +Compute the exponential map on the [`Hyperrectangle`](@ref) manifold `M` from `p` in direction +`X`, which in this case is just +````math +\exp_p X = p + X. +```` +""" +Base.exp(::Hyperrectangle, p, X) = p + X +Base.exp(::Hyperrectangle, p, X, t::Number) = p .+ t .* X + +exp!(::Hyperrectangle, q, p, X) = (q .= p .+ X) +exp!(::Hyperrectangle, q, p, X, t::Number) = (q .= p .+ t .* X) + +function get_coordinates_orthonormal(::Hyperrectangle, p, X, ::RealNumbers) + return vec(X) +end + +function get_coordinates_orthonormal!(::Hyperrectangle, c, p, X, ::RealNumbers) + copyto!(c, vec(X)) + return c +end + +function get_vector_orthonormal(::Hyperrectangle{<:AbstractVector}, ::Any, c, ::RealNumbers) + # this method is defined just to skip a reshape + return c +end + +function get_vector_orthonormal!( + ::Hyperrectangle{<:AbstractVector}, + Y, + ::Any, + c, + ::RealNumbers, +) + # this method is defined just to skip a reshape + copyto!(Y, c) + return Y +end +function get_vector_orthonormal!(M::Hyperrectangle, Y, ::Any, c, ::RealNumbers) + S = representation_size(M) + copyto!(Y, reshape(c, S)) + return Y +end + +@doc raw""" + injectivity_radius(M::Hyperrectangle, p) + +Return the injectivity radius on the [`Hyperrectangle`](@ref) `M` at point `p`, which is +the distance to the nearest boundary the point is not on. +""" +function injectivity_radius(M::Hyperrectangle, p) + ir = Inf + for i in eachindex(M.lb, p) + dist_ub = M.ub[i] - p[i] + if dist_ub > 0 + ir = min(ir, dist_ub) + end + dist_lb = p[i] - M.lb[i] + if dist_lb > 0 + ir = min(ir, dist_lb) + end + end + return ir +end +injectivity_radius(M::Hyperrectangle) = 0.0 +injectivity_radius(M::Hyperrectangle, ::AbstractRetractionMethod) = 0.0 +injectivity_radius(M::Hyperrectangle, p, ::ProjectionRetraction) = injectivity_radius(M, p) +injectivity_radius(M::Hyperrectangle, p, ::ExponentialRetraction) = injectivity_radius(M, p) + +@doc raw""" + inner(M::Hyperrectangle, p, X, Y) + +Compute the inner product on the [`Hyperrectangle`](@ref) `M`, which is just +the inner product on the real-valued vector space of arrays (or tensors) +of size ``n_1 × n_2 × … × n_i``, i.e. + +````math +g_p(X,Y) = \sum_{k ∈ I} X_{k} Y_{k}, +```` + +where ``I`` is the set of vectors ``k ∈ ℕ^i``, such that for all + +``i ≤ j ≤ i`` it holds ``1 ≤ k_j ≤ n_j``. + +For the special case of ``i ≤ 2``, i.e. matrices and vectors, this simplifies to + +````math +g_p(X,Y) = \operatorname{tr}(X^{\mathrm{T}}Y), +```` + +where ``⋅^{\mathrm{T}}`` denotes transposition. +""" +@inline inner(::Hyperrectangle, p, X, Y) = dot(X, Y) +""" + is_flat(::Hyperrectangle) + +Return true. [`Hyperrectangle`](@ref) is a flat manifold. +""" +is_flat(M::Hyperrectangle) = true + +@doc raw""" + log(M::Hyperrectangle, p, q) + +Compute the logarithmic map on the [`Hyperrectangle`](@ref) `M` from `p` to `q`, +which in this case is just +````math +\log_p q = q-p. +```` +""" +Base.log(::Hyperrectangle, ::Any...) +Base.log(::Hyperrectangle, p, q) = q .- p + +log!(::Hyperrectangle, X, p, q) = (X .= q .- p) + +_product_of_dimensions(M::Hyperrectangle) = prod(size(M.lb)) + +""" + manifold_dimension(M::Hyperrectangle) + +Return the manifold dimension of the [`Hyperrectangle`](@ref) `M`, i.e. the product of all array +dimensions. +""" +function manifold_dimension(M::Hyperrectangle) + return _product_of_dimensions(M) +end + +""" + manifold_volume(::Hyperrectangle) + +Return volume of the [`Hyperrectangle`](@ref) manifold, i.e. infinity. +""" +function manifold_volume(M::Hyperrectangle) + vol = one(eltype(M.lb)) + for i in eachindex(M.lb, M.ub) + vol *= M.ub[i] - M.lb[i] + end + return vol +end + +# +# When Statistics / Statsbase.mean! is consistent with mean, we can pass this on to them as well +function Statistics.mean!( + ::Hyperrectangle, + y, + x::AbstractVector, + ::EfficientEstimator; + kwargs..., +) + n = length(x) + copyto!(y, first(x)) + @inbounds for j in 2:n + y .+= x[j] + end + y ./= n + return y +end +function Statistics.mean!( + ::Hyperrectangle, + y, + x::AbstractVector, + w::AbstractWeights, + ::EfficientEstimator; + kwargs..., +) + n = length(x) + if length(w) != n + throw( + DimensionMismatch( + "The number of weights ($(length(w))) does not match the number of points for the mean ($(n)).", + ), + ) + end + copyto!(y, first(x)) + y .*= first(w) + @inbounds for j in 2:n + iszero(w[j]) && continue + y .+= w[j] .* x[j] + end + y ./= sum(w) + return y +end + +mid_point(::Hyperrectangle, p1, p2) = (p1 .+ p2) ./ 2 + +function mid_point!(::Hyperrectangle, q, p1, p2) + q .= (p1 .+ p2) ./ 2 + return q +end + +@doc raw""" + norm(M::Hyperrectangle, p, X) + +Compute the norm of a tangent vector `X` at `p` on the [`Hyperrectangle`](@ref) +`M`, i.e. since every tangent space can be identified with `M` itself +in this case, just the (Frobenius) norm of `X`. +""" +LinearAlgebra.norm(::Hyperrectangle, ::Any, X) = norm(X) + +""" + parallel_transport_direction(M::Hyperrectangle, p, X, d) + +the parallel transport on [`Hyperrectangle`](@ref) is the identiy, i.e. returns `X`. +""" +parallel_transport_direction(::Hyperrectangle, ::Any, X, ::Any) = X +parallel_transport_direction!(::Hyperrectangle, Y, ::Any, X, ::Any) = copyto!(Y, X) + +""" + parallel_transport_to(M::Hyperrectangle, p, X, q) + +the parallel transport on [`Hyperrectangle`](@ref) is the identiy, i.e. returns `X`. +""" +parallel_transport_to(::Hyperrectangle, ::Any, X, ::Any) = X +parallel_transport_to!(::Hyperrectangle, Y, ::Any, X, ::Any) = copyto!(Y, X) + +@doc raw""" + project(M::Hyperrectangle, p) + +Project an arbitrary point `p` onto the [`Hyperrectangle`](@ref) manifold `M`, which +is of course just the identity map. +""" +project(::Hyperrectangle, ::Any) + +function project!(M::Hyperrectangle, q, p) + copyto!(q, p) + for i in eachindex(M.lb, q) + q[i] = clamp(q[i], M.lb[i], M.ub[i]) + end + return q +end + +""" + project(M::Hyperrectangle, p, X) + +Project an arbitrary vector `X` into the tangent space of a point `p` on the +[`Hyperrectangle`](@ref) `M`, which is just the identity, since any tangent +space of `M` can be identified with all of `M`. +""" +project(::Hyperrectangle, ::Any, ::Any) + +function project!(M::Hyperrectangle, Y, p, X) + copyto!(Y, X) + for i in eachindex(M.lb, Y) + if Y[i] >= 0 + Y[i] = min(Y[i], M.ub[i] - p[i]) + else + Y[i] = max(Y[i], M.lb[i] - p[i]) + end + end + return Y +end + +function Random.rand!( + rng::AbstractRNG, + M::Hyperrectangle, + pX; + σ=one(eltype(pX)), + vector_at=nothing, +) + if vector_at === nothing + pX .= M.lb .+ rand(rng, eltype(M.lb), size(M.lb)) .* (M.ub .- M.lb) + else + pX .= randn(rng, eltype(pX), size(pX)) .* σ + project!(M, pX, vector_at, pX) + end + return pX +end + +""" + representation_size(M::Hyperrectangle) + +Return the array dimensions required to represent an element on the +[`Hyperrectangle`](@ref) `M`, i.e. the vector of all array dimensions. +""" +representation_size(M::Hyperrectangle) = size(M.lb) + +function retract_project!(M::Hyperrectangle, r, q, Y, t::Number) + r .= q .+ t .* Y + project!(M, r, r) + return r +end + +@doc raw""" + riemann_tensor(M::Hyperrectangle, p, X, Y, Z) + +Compute the Riemann tensor ``R(X,Y)Z`` at point `p` on [`Hyperrectangle`](@ref) manifold `M`. +Its value is always the zero tangent vector. +""" +riemann_tensor(M::Hyperrectangle, p, X, Y, Z) + +function riemann_tensor!(::Hyperrectangle, Xresult, p, X, Y, Z) + return fill!(Xresult, 0) +end + +@doc raw""" + sectional_curvature(::Hyperrectangle, p, X, Y) + +Sectional curvature of [`Hyperrectangle`](@ref) manifold `M` is 0. +""" +function sectional_curvature(::Hyperrectangle, p, X, Y) + return 0.0 +end + +@doc raw""" + sectional_curvature_max(::Hyperrectangle) + +Sectional curvature of [`Hyperrectangle`](@ref) manifold `M` is 0. +""" +function sectional_curvature_max(::Hyperrectangle) + return 0.0 +end + +@doc raw""" + sectional_curvature_min(M::Hyperrectangle) + +Sectional curvature of [`Hyperrectangle`](@ref) manifold `M` is 0. +""" +function sectional_curvature_min(::Hyperrectangle) + return 0.0 +end + +function Base.show(io::IO, M::Hyperrectangle) + return print(io, "Hyperrectangle($(M.lb), $(M.ub))") +end + +function vector_transport_direction( + M::Hyperrectangle, + p, + X, + ::Any, + ::AbstractVectorTransportMethod=default_vector_transport_method(M, typeof(p)), + ::AbstractRetractionMethod=default_retraction_method(M, typeof(p)), +) + return X +end +function vector_transport_direction!( + M::Hyperrectangle, + Y, + p, + X, + ::Any, + ::AbstractVectorTransportMethod=default_vector_transport_method(M, typeof(p)), + ::AbstractRetractionMethod=default_retraction_method(M, typeof(p)), +) + return copyto!(Y, X) +end +""" + vector_transport_to(M::Hyperrectangle, p, X, q, ::AbstractVectorTransportMethod) + +Transport the vector `X` from the tangent space at `p` to the tangent space at `q` +on the [`Hyperrectangle`](@ref) `M`, which simplifies to the identity. +""" +vector_transport_to(::Hyperrectangle, ::Any, ::Any, ::Any, ::AbstractVectorTransportMethod) +function vector_transport_to( + M::Hyperrectangle, + p, + X, + ::Any, + ::AbstractVectorTransportMethod=default_vector_transport_method(M, typeof(p)), + ::AbstractRetractionMethod=default_retraction_method(M, typeof(p)), +) + return X +end + +function vector_transport_to!( + M::Hyperrectangle, + Y, + p, + X, + ::Any, + ::AbstractVectorTransportMethod=default_vector_transport_method(M, typeof(p)), + ::AbstractRetractionMethod=default_retraction_method(M, typeof(p)), +) + return copyto!(Y, X) +end + +Statistics.var(::Hyperrectangle, x::AbstractVector; kwargs...) = sum(var(x; kwargs...)) + +@doc raw""" + volume_density(M::Hyperrectangle, p, X) + +Return volume density function of [`Hyperrectangle`](@ref) manifold `M`, i.e. 1. +""" +function volume_density(::Hyperrectangle, p, X) + return one(eltype(X)) +end + +@doc raw""" + Y = Weingarten(M::Hyperrectangle, p, X, V) + Weingarten!(M::Hyperrectangle, Y, p, X, V) + +Compute the Weingarten map ``\mathcal W_p`` at `p` on the [`Hyperrectangle`](@ref) `M` with +respect to the tangent vector ``X \in T_p\mathcal M`` and the normal vector ``V \in N_p\mathcal M``. + +Since this a flat space by itself, the result is always the zero tangent vector. +""" +Weingarten(::Hyperrectangle, p, X, V) + +Weingarten!(::Hyperrectangle, Y, p, X, V) = fill!(Y, 0) + +""" + zero_vector(M::Hyperrectangle, p) + +Return the zero vector in the tangent space of `p` on the [`Hyperrectangle`](@ref) +`M`, which here is just a zero filled array the same size as `p`. +""" +zero_vector(::Hyperrectangle, ::Any...) + +zero_vector!(::Hyperrectangle, X, ::Any) = fill!(X, 0) diff --git a/test/manifolds/circle.jl b/test/manifolds/circle.jl index 8296554eb0..9c3b96765f 100644 --- a/test/manifolds/circle.jl +++ b/test/manifolds/circle.jl @@ -156,6 +156,7 @@ using Manifolds: TFVector, CoTFVector is_mutating=false, test_rand_point=true, test_rand_tvector=true, + rand_tvector_atol_multiplier=2.0, ) ptsS = map(p -> (@SArray fill(p)), pts) test_manifold( @@ -175,6 +176,7 @@ using Manifolds: TFVector, CoTFVector basis_types_to_from=basis_types_real, test_rand_point=true, test_rand_tvector=true, + rand_tvector_atol_multiplier=2.0, ) end end diff --git a/test/manifolds/hyperrectangle.jl b/test/manifolds/hyperrectangle.jl new file mode 100644 index 0000000000..43a8868e3b --- /dev/null +++ b/test/manifolds/hyperrectangle.jl @@ -0,0 +1,136 @@ +include("../header.jl") + +@testset "Hyperrectangle" begin + M = Hyperrectangle([-1.0, 2.0, -3.0], [1.0, 4.0, 9.0]) + @test repr(M) == "Hyperrectangle([-1.0, 2.0, -3.0], [1.0, 4.0, 9.0])" + + @test_throws ArgumentError Hyperrectangle([1.0, 2.0], [-1.0, 5.0]) + + @test is_flat(M) + p = zeros(3) + @test project!(M, p, p) == p + @test embed!(M, p, p) == p + X = zeros(3) + X[1] = 1.0 + Y = similar(X) + project!(M, Y, p, X) + @test Y == X + @test embed(M, p, X) == X + + @test_throws DomainError is_point(M, [1.0im, 0.0, 0.0]; error=:error) + @test_throws DomainError is_point(M, [10.0, 3.0, 0.0]; error=:error) + @test_throws DomainError is_vector(M, [1.0, 2.0, 0.0], [1.0im, 0.0, 0.0]; error=:error) + @test_throws DomainError is_vector(M, [1], [1.0, 1.0, 0.0]; error=:error) + @test_throws DomainError is_vector(M, [0.0, 0.0, 0.0], [1.0]; error=:error) + @test_throws DomainError is_vector(M, [0.0, 0.0, 0.0], [1.0, 0.0, 1.0im]; error=:error) + + @testset "projections" begin + @test project(M, [4.0, -2.0, 3.0]) ≈ [1.0, 2.0, 3.0] + @test project(M, [1.0, 2.0, 3.0], [2.0, 0.5, -10.0]) ≈ [0.0, 0.5, -6.0] + end + + @testset "injectivity_radius" begin + @test injectivity_radius(M) == 0.0 + @test injectivity_radius(M, ExponentialRetraction()) == 0.0 + @test injectivity_radius(M, [0.0, 2.5, 1.0]) == 0.5 + @test injectivity_radius(M, [0.0, 2.5, 1.0], ExponentialRetraction()) == 0.5 + @test injectivity_radius(M, [0.0, 2.5, 1.0], ProjectionRetraction()) == 0.5 + @test injectivity_radius(M, [0.0, 2.0, 1.0]) == 1.0 + @test injectivity_radius(M, [0.5, 4.0, 1.0]) == 0.5 + + M2 = Hyperrectangle([-1.0, 2.0, -3.0], [1.0, 2.0, 9.0]) + @test injectivity_radius(M2, [0.0, 2.0, 1.0]) == 1.0 + end + + basis_types = (DefaultOrthonormalBasis(),) + pts = [[1.0, 2.0, 0.0], [0.0, 3.0, 0.0], [0.0, 3.5, 1.0]] + + test_manifold( + M, + pts, + test_injectivity_radius=false, + parallel_transport=true, + test_project_point=true, + test_project_tangent=true, + test_default_vector_transport=true, + vector_transport_methods=[ParallelTransport()], + retraction_methods=[ProjectionRetraction()], + test_mutating_rand=true, + basis_types_vecs=basis_types, + basis_types_to_from=basis_types, + test_inplace=true, + test_rand_point=true, + test_rand_tvector=true, + ) + @testset "Array Hyperrectangle" begin + MA = Hyperrectangle([-1.0 2.0 3.0; 2.0 5.0 10.0], [10.0 12.0 13.0; 12.0 15.0 20.0]) + pts_a = [ + [2.0 2.0 3.0; 2.0 6.0 10.0], + [-1.0 5.0 3.0; 4.0 5.0 10.0], + [-1.0 2.0 5.0; 6.0 6.0 15.0], + ] + test_manifold( + MA, + pts_a, + test_injectivity_radius=false, + test_project_point=true, + test_project_tangent=true, + test_default_vector_transport=true, + vector_transport_methods=[ParallelTransport()], + test_mutating_rand=true, + basis_types_vecs=basis_types, + basis_types_to_from=basis_types, + test_inplace=true, + test_rand_point=true, + test_rand_tvector=true, + ) + end + + @testset "Hyperrectangle(1)" begin + M = Hyperrectangle([0.0], [5.0]) + @test distance(M, 2.0, 4.0) == 2.0 + end + + @testset "errors" begin + M = Hyperrectangle([-1.0, 2.0, -3.0, 0.0], [1.0, 4.0, 9.0, 3.0]) + @test_throws DimensionMismatch distance(M, [1, 2, 3, 4], [1 2; 3 4]) + end + + @testset "Weingarten & Hessian" begin + M = Hyperrectangle([-1.0, -1.0], [10.0, 10.0]) + p = [1.0, 2.0] + G = [3.0, 4.0] + H = [5.0, 6.0] + X = [7.0, 8.0] + rH = riemannian_Hessian(M, p, G, H, X) + @test rH == H + end + @testset "Volume" begin + M = Hyperrectangle([-1.0, 2.0, -3.0], [1.0, 4.0, 9.0]) + @test manifold_volume(M) == 48.0 + p = zeros(3) + X = zeros(3) + @test volume_density(M, p, X) == 1.0 + end + + @testset "statistics" begin + @test default_approximation_method(M, mean) === EfficientEstimator() + @test mean(M, pts) ≈ [1 / 3, 17 / 6, 1 / 3] + @test mean(M, pts, pweights(ones(3) / 3)) ≈ [1 / 3, 17 / 6, 1 / 3] + @test_throws DimensionMismatch mean(M, pts, pweights(ones(4) / 4)) + @test var(M, pts) ≈ 1.25 + end + + @testset "Euclidean metric tests" begin + @test riemann_tensor( + M, + pts[1], + pts[2] - pts[1], + pts[3] - pts[1], + (pts[3] - pts[1]) / 2, + ) == [0.0, 0.0, 0.0] + @test sectional_curvature(M, p, [1.0, 0.0], [0.0, 1.0]) == 0.0 + @test sectional_curvature_max(M) == 0.0 + @test sectional_curvature_min(M) == 0.0 + end +end diff --git a/test/manifolds/sphere.jl b/test/manifolds/sphere.jl index e81c340079..643f99c968 100644 --- a/test/manifolds/sphere.jl +++ b/test/manifolds/sphere.jl @@ -308,4 +308,23 @@ using ManifoldsBase: TFVector p = [1.0, 0.0, 0.0] @test local_metric(M, p, DefaultOrthonormalBasis()) == Diagonal([1.0, 1.0]) end + + @testset "sectional curvature" begin + M = Sphere(2; parameter=:field) + K = Manifolds.sectional_curvature_matrix( + M, + [1.0, 0.0, 0.0], + DefaultOrthonormalBasis(), + ) + @test isapprox(K, [0 1; 1 0]) + @test isapprox( + Manifolds.estimated_sectional_curvature_matrix( + M, + [1.0, 0.0, 0.0], + DefaultOrthonormalBasis(), + ), + [0.0 1.0; 1.0 0.0], + atol=0.15, + ) + end end diff --git a/test/manifolds/symmetric_positive_definite.jl b/test/manifolds/symmetric_positive_definite.jl index b01662de22..6bfc0c282f 100644 --- a/test/manifolds/symmetric_positive_definite.jl +++ b/test/manifolds/symmetric_positive_definite.jl @@ -313,4 +313,10 @@ include("../header.jl") @test typeof(get_embedding(M)) === Euclidean{Tuple{Int,Int},ℝ} @test repr(M) == "SymmetricPositiveDefinite(3; parameter=:field)" end + + @testset "Curvature" begin + @test sectional_curvature_min(SymmetricPositiveDefinite(1)) == 0.0 + @test sectional_curvature_min(SymmetricPositiveDefinite(3)) == -0.25 + @test sectional_curvature_max(SymmetricPositiveDefinite(3)) == 0.0 + end end diff --git a/test/runtests.jl b/test/runtests.jl index 1764273c24..f6bd62c799 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -151,6 +151,7 @@ end include_test("manifolds/grassmann.jl") include_test("manifolds/hamiltonian.jl") include_test("manifolds/hyperbolic.jl") + include_test("manifolds/hyperrectangle.jl") include_test("manifolds/lorentz.jl") include_test("manifolds/multinomial_doubly_stochastic.jl") include_test("manifolds/multinomial_symmetric.jl")