From deede80e418992882ae51ea9e1c0878cd28cb86d Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Mon, 25 Dec 2023 08:25:24 +0100 Subject: [PATCH 1/5] Unify use of an `atol` when comparing to zero (#692) * Unify atols in isapproxes. * fix dispatch route for direction vector transport in Circle. * Make the atol more generic/tolerant * update the literature entry for the Manifolds.jl paper --------- Co-authored-by: Mateusz Baran --- NEWS.md | 9 + Project.toml | 6 +- README.md | 18 +- docs/src/references.bib | 16 +- src/Manifolds.jl | 32 +- src/deprecated.jl | 5 + src/groups/addition_operation.jl | 10 +- src/groups/group_action.jl | 6 +- src/groups/group_operation_action.jl | 2 +- src/groups/special_linear.jl | 10 +- src/manifolds/CenteredMatrices.jl | 19 +- src/manifolds/CholeskySpace.jl | 19 +- src/manifolds/Circle.jl | 17 +- src/manifolds/Elliptope.jl | 10 +- src/manifolds/EmbeddedTorus.jl | 4 +- src/manifolds/EssentialManifold.jl | 10 +- src/manifolds/Euclidean.jl | 63 +++- src/manifolds/FixedRankMatrices.jl | 12 +- src/manifolds/GeneralUnitaryMatrices.jl | 2 +- src/manifolds/GeneralizedGrassmann.jl | 2 +- src/manifolds/Grassmann.jl | 2 +- src/manifolds/Hyperbolic.jl | 2 +- src/manifolds/HyperbolicHyperboloid.jl | 10 +- src/manifolds/KendallsPreShapeSpace.jl | 19 +- src/manifolds/MultinomialDoublyStochastic.jl | 19 +- src/manifolds/ProbabilitySimplex.jl | 12 +- src/manifolds/ProjectiveSpace.jl | 12 +- src/manifolds/Rotations.jl | 12 +- src/manifolds/SPDFixedDeterminant.jl | 10 +- src/manifolds/Spectrahedron.jl | 10 +- src/manifolds/Sphere.jl | 12 +- src/manifolds/SphereSymmetricMatrices.jl | 19 +- src/manifolds/Symmetric.jl | 8 +- src/manifolds/SymmetricPositiveDefinite.jl | 4 +- .../SymmetricPositiveSemidefiniteFixedRank.jl | 12 +- src/manifolds/Symplectic.jl | 19 +- src/manifolds/SymplecticStiefel.jl | 21 +- src/statistics.jl | 296 +++++++----------- test/manifolds/circle.jl | 6 + test/manifolds/embedded_torus.jl | 2 +- test/manifolds/rotations.jl | 29 ++ test/statistics.jl | 26 +- test/test_deprecated.jl | 34 +- 43 files changed, 558 insertions(+), 310 deletions(-) diff --git a/NEWS.md b/NEWS.md index e54201ab67..89c862fe12 100644 --- a/NEWS.md +++ b/NEWS.md @@ -5,6 +5,15 @@ 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.9] – 2023-12-24 + +### Fixed + +* introduced a nonzero `atol` for all point and vector checks that compre to zero. + This makes those checks a bit more relaxed by default and resolves [#630](https://github.com/JuliaManifolds/Manifolds.jl/issues/630). +* `default_estimation_method(M, f)` is deprecated, use `default_approximation_method(M, f)` for your specific method `f` on the manifold `M`. +* `AbstractEstimationMethod` is deprecated, use `AbstractApproximationMethod` instead. + ## [0.9.8] - 2023-11-17 ### Fixed diff --git a/Project.toml b/Project.toml index 490fd5c2bf..54bb294f41 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.8" +version = "0.9.9" [deps] Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" @@ -42,7 +42,7 @@ ManifoldsRecipesBaseExt = ["Colors", "RecipesBase"] ManifoldsTestExt = "Test" [compat] -BoundaryValueDiffEq = "4, 5" +BoundaryValueDiffEq = "4, 5.6.1" Colors = "0.12" Distributions = "0.22.6, 0.23, 0.24, 0.25" Einsum = "0.4" @@ -51,7 +51,7 @@ HybridArrays = "0.4" Kronecker = "0.4, 0.5" LinearAlgebra = "1.6" ManifoldDiff = "0.3.7" -ManifoldsBase = "0.15.0" +ManifoldsBase = "0.15.6" Markdown = "1.6" MatrixEquations = "2.2" OrdinaryDiffEq = "6.31" diff --git a/README.md b/README.md index afd446fba5..8cebb371b0 100644 --- a/README.md +++ b/README.md @@ -53,17 +53,19 @@ If you have any questions regarding the Manifolds.jl ecosystem feel free to reac ## Citation -If you use `Manifolds.jl` in your work, please cite the following +If you use `Manifolds.jl` in your work, please cite the following open access article ```biblatex @article{AxenBaranBergmannRzecki:2023, - AUTHOR = {Seth D. Axen and Mateusz Baran and Ronny Bergmann and Krzysztof Rzecki}, - EPRINT = {2021.08777}, - EPRINTTYPE = {arXiv}, - JOURNAL = {AMS Transactions on Mathematical Software}, - NOTE = {accepted for publication}, - TITLE = {Manifolds.jl: An Extensible {J}ulia Framework for Data Analysis on Manifolds}, - YEAR = {2023} + AUTHOR = {Axen, Seth D. and Baran, Mateusz and Bergmann, Ronny and Rzecki, Krzysztof}, + ARTICLENO = {33}, + DOI = {10.1145/3618296}, + JOURNAL = {ACM Transactions on Mathematical Software}, + MONTH = {dec}, + NUMBER = {4}, + TITLE = {Manifolds.Jl: An Extensible Julia Framework for Data Analysis on Manifolds}, + VOLUME = {49}, + YEAR = {2023} } ``` diff --git a/docs/src/references.bib b/docs/src/references.bib index ecd9632754..0b43a860f6 100644 --- a/docs/src/references.bib +++ b/docs/src/references.bib @@ -89,13 +89,15 @@ @book{AyJostLeSchwachhoefer:2017 PUBLISHER = {Springer Cham} } @article{AxenBaranBergmannRzecki:2023, - AUTHOR = {Seth D. Axen and Mateusz Baran and Ronny Bergmann and Krzysztof Rzecki}, - EPRINT = {2021.08777}, - EPRINTTYPE = {arXiv}, - JOURNAL = {AMS Transactions on Mathematical Software}, - NOTE = {accepted for publication}, - TITLE = {Manifolds.jl: An Extensible {J}ulia Framework for Data Analysis on Manifolds}, - YEAR = {2023} + AUTHOR = {Axen, Seth D. and Baran, Mateusz and Bergmann, Ronny and Rzecki, Krzysztof}, + ARTICLENO = {33}, + DOI = {10.1145/3618296}, + JOURNAL = {ACM Transactions on Mathematical Software}, + MONTH = {dec}, + NUMBER = {4}, + TITLE = {Manifolds.Jl: An Extensible Julia Framework for Data Analysis on Manifolds}, + VOLUME = {49}, + YEAR = {2023} } # # B diff --git a/src/Manifolds.jl b/src/Manifolds.jl index 3ca4e4c428..89c466e728 100644 --- a/src/Manifolds.jl +++ b/src/Manifolds.jl @@ -48,6 +48,7 @@ import ManifoldsBase: check_vector, copy, copyto!, + default_approximation_method, default_inverse_retraction_method, default_retraction_method, default_vector_transport_method, @@ -149,19 +150,14 @@ import ManifoldsBase: submanifold_component, submanifold_components, vector_space_dimension, - vector_transport_along, # just specified in Euclidean - the next 5 as well - vector_transport_along_diff, - vector_transport_along_project, + vector_transport_along, # just specified in Euclidean - the next 5 as well vector_transport_along!, - vector_transport_along_diff!, - vector_transport_along_project!, + vector_transport_along_diff!, # For consistency these are imported, but for now not + vector_transport_along_project!, # overwritten with new definitons. vector_transport_direction, - vector_transport_direction_diff, vector_transport_direction!, vector_transport_direction_diff!, vector_transport_to, - vector_transport_to_diff, - vector_transport_to_project, vector_transport_to!, vector_transport_to_diff!, vector_transport_to_project!, # some overwrite layer 2 @@ -186,6 +182,9 @@ import ManifoldDiff: riemannian_Hessian, riemannian_Hessian! +import Statistics: mean, mean!, median, median!, cov, var +import StatsBase: mean_and_var + using Base.Iterators: repeated using Distributions using Einsum: @einsum @@ -198,6 +197,7 @@ using ManifoldsBase: ℝ, ℂ, ℍ, + AbstractApproximationMethod, AbstractBasis, AbstractDecoratorManifold, AbstractInverseRetractionMethod, @@ -225,6 +225,7 @@ using ManifoldsBase: CotangentSpaceType, CoTFVector, CoTVector, + CyclicProximalPointEstimation, DefaultBasis, DefaultOrthogonalBasis, DefaultOrthonormalBasis, @@ -232,13 +233,18 @@ using ManifoldsBase: DiagonalizingBasisData, DiagonalizingOrthonormalBasis, DifferentiatedRetractionVectorTransport, + EfficientEstimator, EmbeddedManifold, EmptyTrait, EuclideanMetric, ExponentialRetraction, + ExtrinsicEstimation, Fiber, FiberType, FVector, + GeodesicInterpolation, + GeodesicInterpolationWithinRadius, + GradientDescentEstimation, InverseProductRetraction, IsIsometricEmbeddedManifold, IsEmbeddedManifold, @@ -295,6 +301,7 @@ using ManifoldsBase: VectorSpaceFiber, VectorSpaceType, VeeOrthogonalBasis, + WeiszfeldEstimation, @invoke_maker, _euclidean_basis_vector, combine_allocation_promotion_functions, @@ -607,8 +614,9 @@ include("deprecated.jl") export test_manifold export test_group, test_action -# +# Abstract main types export CoTVector, AbstractManifold, AbstractManifoldPoint, TVector +# Manifolds export AbstractSphere, AbstractProjectiveSpace export Euclidean, ArrayProjectiveSpace, @@ -748,9 +756,10 @@ export AbstractInverseRetractionMethod, ShootingInverseRetraction, SoftmaxInverseRetraction # Estimation methods for median and mean -export AbstractEstimationMethod, +export AbstractApproximationMethod, GradientDescentEstimation, CyclicProximalPointEstimation, + EfficientEstimator, GeodesicInterpolation, GeodesicInterpolationWithinRadius, ExtrinsicEstimation @@ -788,9 +797,10 @@ export ×, convert, complex_dot, decorated_manifold, - default_vector_transport_method, + default_approximation_method, default_inverse_retraction_method, default_retraction_method, + default_vector_transport_method, det_local_metric, differential_canonical_project, differential_canonical_project!, diff --git a/src/deprecated.jl b/src/deprecated.jl index 8b13789179..03601238b0 100644 --- a/src/deprecated.jl +++ b/src/deprecated.jl @@ -1 +1,6 @@ +@deprecate default_estimation_method(M::AbstractManifold, f) default_approximation_method( + M, + f, +) +@deprecate ExtrinsicEstimation() ExtrinsicEstimation(EfficientEstimator()) diff --git a/src/groups/addition_operation.jl b/src/groups/addition_operation.jl index d4a8920e7c..885974cf53 100644 --- a/src/groups/addition_operation.jl +++ b/src/groups/addition_operation.jl @@ -83,8 +83,14 @@ function inv_diff!(::AdditionGroupTrait, G::AbstractDecoratorManifold, Y, p, X) return Y end -function is_identity(::AdditionGroupTrait, G::AbstractDecoratorManifold, q; kwargs...) - return isapprox(G, q, zero(q); kwargs...) +function is_identity( + ::AdditionGroupTrait, + G::AbstractDecoratorManifold, + q::T; + atol::Real=sqrt(prod(representation_size(G))) * eps(real(float(number_eltype(T)))), + kwargs..., +) where {T} + return isapprox(G, q, zero(q); atol=atol, kwargs...) end # resolve ambiguities function is_identity( diff --git a/src/groups/group_action.jl b/src/groups/group_action.jl index 782ce45ecb..d31dd18786 100644 --- a/src/groups/group_action.jl +++ b/src/groups/group_action.jl @@ -116,7 +116,7 @@ where element `a` is acting on `p`, with respect to the group element. Let ``\mathcal G`` be the group acting on manifold ``\mathcal M`` by the action `A`. The action is of element ``g ∈ \mathcal G`` on a point ``p ∈ \mathcal M``. -The differential transforms vector `X` from the tangent space at `a ∈ \mathcal G`, +The differential transforms vector `X` from the tangent space at `a ∈ \mathcal G`, ``X ∈ T_a \mathcal G`` into a tangent space of the manifold ``\mathcal M``. When action on element `p` is written as ``\mathrm{d}τ^p``, with the specified left or right convention, the differential transforms vectors @@ -193,7 +193,7 @@ end A::AbstractGroupAction, pts, p, - mean_method::AbstractEstimationMethod = GradientDescentEstimation(), + mean_method::AbstractApproximationMethod = GradientDescentEstimation(), ) Calculate an action element ``a`` of action `A` that is the mean element of the orbit of `p` @@ -210,7 +210,7 @@ function center_of_orbit( A::AbstractGroupAction, pts::AbstractVector, q, - mean_method::AbstractEstimationMethod=GradientDescentEstimation(), + mean_method::AbstractApproximationMethod=GradientDescentEstimation(), ) alignments = map(p -> optimal_alignment(A, q, p), pts) return mean(base_group(A), alignments, mean_method) diff --git a/src/groups/group_operation_action.jl b/src/groups/group_operation_action.jl index 2a7d06557d..1b9098ed15 100644 --- a/src/groups/group_operation_action.jl +++ b/src/groups/group_operation_action.jl @@ -199,7 +199,7 @@ function center_of_orbit( A::GroupOperationAction, pts::AbstractVector, p, - mean_method::AbstractEstimationMethod, + mean_method::AbstractApproximationMethod, ) μ = mean(A.group, pts, mean_method) return inverse_apply(switch_direction_and_side(A), p, μ) diff --git a/src/groups/special_linear.jl b/src/groups/special_linear.jl index 7ffb5fde1e..e8f8eb05a9 100644 --- a/src/groups/special_linear.jl +++ b/src/groups/special_linear.jl @@ -49,9 +49,15 @@ function check_point(G::SpecialLinear, p; kwargs...) return nothing end -function check_vector(G::SpecialLinear, p, X; kwargs...) +function check_vector( + G::SpecialLinear, + p, + X::T; + atol::Real=sqrt(prod(representation_size(G))) * eps(real(float(number_eltype(T)))), + kwargs..., +) where {T} trX = tr(inverse_translate_diff(G, p, p, X, LeftForwardAction())) - if !isapprox(trX, 0; kwargs...) + if !isapprox(trX, 0; atol=atol, kwargs...) return DomainError( trX, "The matrix $(X) does not lie in the tangent space of $(G) at $(p), since " * diff --git a/src/manifolds/CenteredMatrices.jl b/src/manifolds/CenteredMatrices.jl index 0856c49a52..e5f801797a 100644 --- a/src/manifolds/CenteredMatrices.jl +++ b/src/manifolds/CenteredMatrices.jl @@ -35,9 +35,14 @@ zero. The tolerance for the column sums of `p` can be set using `kwargs...`. """ -function check_point(M::CenteredMatrices, p; kwargs...) +function check_point( + M::CenteredMatrices, + p::T; + atol::Real=sqrt(prod(representation_size(M))) * eps(real(float(number_eltype(T)))), + kwargs..., +) where {T} m, n = get_parameter(M.size) - if !isapprox(sum(p, dims=1), zeros(1, n); kwargs...) + if !isapprox(sum(p, dims=1), zeros(1, n); atol=atol, kwargs...) return DomainError( p, string( @@ -56,9 +61,15 @@ Check whether `X` is a tangent vector to manifold point `p` on the sum to zero and its values are from the correct [`AbstractNumbers`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/types.html#number-system). The tolerance for the column sums of `p` and `X` can be set using `kwargs...`. """ -function check_vector(M::CenteredMatrices, p, X; kwargs...) +function check_vector( + M::CenteredMatrices, + p, + X::T; + atol::Real=sqrt(prod(representation_size(M))) * eps(real(float(number_eltype(T)))), + kwargs..., +) where {T} m, n = get_parameter(M.size) - if !isapprox(sum(X, dims=1), zeros(1, n); kwargs...) + if !isapprox(sum(X, dims=1), zeros(1, n); atol=atol, kwargs...) return DomainError( X, "The vector $(X) is not a tangent vector to $(p) on $(M), since its columns do not sum to zero.", diff --git a/src/manifolds/CholeskySpace.jl b/src/manifolds/CholeskySpace.jl index 5127235a07..2d2226c4b4 100644 --- a/src/manifolds/CholeskySpace.jl +++ b/src/manifolds/CholeskySpace.jl @@ -28,10 +28,15 @@ it's size fits the manifold, it is a lower triangular matrix and has positive entries on the diagonal. The tolerance for the tests can be set using the `kwargs...`. """ -function check_point(M::CholeskySpace, p; kwargs...) +function check_point( + M::CholeskySpace, + p::T; + atol::Real=sqrt(prod(representation_size(M))) * eps(real(float(number_eltype(T)))), + kwargs..., +) where {T} cks = check_size(M, p) cks === nothing || return cks - if !isapprox(norm(strictlyUpperTriangular(p)), 0.0; kwargs...) + if !isapprox(norm(strictlyUpperTriangular(p)), 0.0; atol=atol, kwargs...) return DomainError( norm(UpperTriangular(p) - Diagonal(p)), "The point $(p) does not lie on $(M), since it strictly upper triangular nonzero entries", @@ -54,8 +59,14 @@ after [`check_point`](@ref)`(M,p)`, `X` has to have the same dimension as `p` and a symmetric matrix. The tolerance for the tests can be set using the `kwargs...`. """ -function check_vector(M::CholeskySpace, p, X; kwargs...) - if !isapprox(norm(strictlyUpperTriangular(X)), 0.0; kwargs...) +function check_vector( + M::CholeskySpace, + p, + X; + atol::Real=sqrt(prod(representation_size(M)) * eps(float(eltype(p)))), + kwargs..., +) + if !isapprox(norm(strictlyUpperTriangular(X)), 0.0; atol=atol, kwargs...) return DomainError( norm(UpperTriangular(X) - Diagonal(X)), "The matrix $(X) is not a tangent vector at $(p) (represented as an element of the Lie algebra) since it is not lower triangular.", diff --git a/src/manifolds/Circle.jl b/src/manifolds/Circle.jl index bacf0320a1..c6b2d7a2f3 100644 --- a/src/manifolds/Circle.jl +++ b/src/manifolds/Circle.jl @@ -78,8 +78,14 @@ check_vector(::Circle{ℝ}, ::Any...; ::Any...) function check_vector(M::Circle{ℝ}, p, X; kwargs...) return nothing end -function check_vector(M::Circle{ℂ}, p, X; kwargs...) - if !isapprox(abs(complex_dot(p, X)), 0.0; kwargs...) +function check_vector( + M::Circle{ℂ}, + p, + X::T; + atol::Real=sqrt(eps(real(float(number_eltype(T))))), + kwargs..., +) where {T} + if !isapprox(abs(complex_dot(p, X)), 0; atol=atol, kwargs...) return DomainError( abs(complex_dot(p, X)), "The value $(X) is not a tangent vector to $(p) on $(M), since it is not orthogonal in the embedding.", @@ -160,7 +166,7 @@ function Base.exp(M::Circle{ℂ}, p::Number, X::Number, t::Number) end exp!(::Circle{ℝ}, q, p, X) = (q .= sym_rem(p + X)) -exp!(::Circle{ℝ}, q, p, X, t::Number) = (q .= sym_rem(p + t * X)) +exp!(::Circle{ℝ}, q, p, X, t::Number) = (q .= sym_rem(p[] + t * X[])) function exp!(M::Circle{ℂ}, q, p, X) θ = norm(M, p, X) q .= cos(θ) * p + usinc(θ) * X @@ -538,6 +544,11 @@ function parallel_transport_to!(M::Circle{ℂ}, Y, p, X, q) return Y end +# dispatch before allocation +function _vector_transport_direction(M::Circle, p, X, d, ::ParallelTransport) + return parallel_transport_to(M, p, X, exp(M, p, d)) +end + """ volume_density(::Circle, p, X) diff --git a/src/manifolds/Elliptope.jl b/src/manifolds/Elliptope.jl index 18b428c3af..e7d7d7432d 100644 --- a/src/manifolds/Elliptope.jl +++ b/src/manifolds/Elliptope.jl @@ -85,10 +85,16 @@ zero diagonal. The tolerance for the base point check and zero diagonal can be set using the `kwargs...`. Note that symmetric of $X$ holds by construction an is not explicitly checked. """ -function check_vector(M::Elliptope, q, Y; kwargs...) +function check_vector( + M::Elliptope, + q, + Y::T; + atol::Real=sqrt(prod(representation_size(M))) * eps(real(float(number_eltype(T)))), + kwargs..., +) where {T} X = q * Y' + Y * q' n = diag(X) - if !all(isapprox.(n, 0.0; kwargs...)) + if !all(isapprox.(n, 0.0; atol=atol, kwargs...)) return DomainError( n, "The vector $(X) is not a tangent to a point on $(M) (represented py $(q) and $(Y), since its diagonal is nonzero.", diff --git a/src/manifolds/EmbeddedTorus.jl b/src/manifolds/EmbeddedTorus.jl index ae06b8cd86..9adb27d2d2 100644 --- a/src/manifolds/EmbeddedTorus.jl +++ b/src/manifolds/EmbeddedTorus.jl @@ -52,9 +52,9 @@ Check whether `X` is a valid vector tangent to `p` on the [`EmbeddedTorus`](@ref The method checks if the vector `X` is orthogonal to the vector normal to the torus, see [`normal_vector`](@ref). Absolute tolerance can be set using `atol`. """ -function check_vector(M::EmbeddedTorus, p, X; atol=eps(eltype(p)), kwargs...) +function check_vector(M::EmbeddedTorus, p, X; atol::Real=eps(float(eltype(p))), kwargs...) dot_nX = dot(normal_vector(M, p), X) - if !isapprox(dot_nX, 0; atol, kwargs...) + if !isapprox(dot_nX, 0; atol=atol, kwargs...) return DomainError(dot_nX, "The vector $(X) is not tangent to $(p) from $(M).") end return nothing diff --git a/src/manifolds/EssentialManifold.jl b/src/manifolds/EssentialManifold.jl index 8d46f2ffca..6d585d4974 100644 --- a/src/manifolds/EssentialManifold.jl +++ b/src/manifolds/EssentialManifold.jl @@ -125,8 +125,14 @@ exp(::EssentialManifold, ::Any...) get_iterator(::EssentialManifold) = Base.OneTo(2) -function _isapprox(M::EssentialManifold, p, q; kwargs...) - return isapprox(distance(M, p, q), 0.0; kwargs...) +function _isapprox( + M::EssentialManifold, + p, + q::T; + atol::Real=eps(real(float(number_eltype(number_eltype(T))))), + kwargs..., +) where {T} + return isapprox(distance(M, p, q), 0.0; atol=atol, kwargs...) end """ diff --git a/src/manifolds/Euclidean.jl b/src/manifolds/Euclidean.jl index c5173c21c0..f2160ef492 100644 --- a/src/manifolds/Euclidean.jl +++ b/src/manifolds/Euclidean.jl @@ -114,6 +114,11 @@ function check_vector(M::Euclidean{N,𝔽}, p, X; kwargs...) where {N,𝔽} return nothing end +default_approximation_method(::Euclidean, ::typeof(mean)) = EfficientEstimator() +function default_approximation_method(::Euclidean, ::typeof(median), ::Type{<:Number}) + return EfficientEstimator() +end + function det_local_metric( ::MetricManifold{𝔽,<:AbstractManifold,EuclideanMetric}, p, @@ -504,23 +509,65 @@ Return volume of the [`Euclidean`](@ref) manifold, i.e. infinity. """ manifold_volume(::Euclidean) = Inf -Statistics.mean(::Euclidean{Tuple{}}, x::AbstractVector{<:Number}; kwargs...) = mean(x) function Statistics.mean( ::Union{Euclidean{TypeParameter{Tuple{}}},Euclidean{Tuple{}}}, - x::AbstractVector{<:Number}; + x::AbstractVector, + ::EfficientEstimator; kwargs..., ) return mean(x) end function Statistics.mean( ::Union{Euclidean{TypeParameter{Tuple{}}},Euclidean{Tuple{}}}, - x::AbstractVector{<:Number}, - w::AbstractWeights; + x::AbstractVector, + w::AbstractWeights, + ::EfficientEstimator; kwargs..., ) return mean(x, w) end -Statistics.mean(::Euclidean, x::AbstractVector; kwargs...) = mean(x) +# +# When Statistics / Statsbase.mean! is consistent with mean, we can pass this on to them as well +function Statistics.mean!( + ::Euclidean, + 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!( + ::Euclidean, + 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 function StatsBase.mean_and_var( ::Union{Euclidean{TypeParameter{Tuple{}}},Euclidean{Tuple{}}}, @@ -543,7 +590,8 @@ end function Statistics.median( ::Union{Euclidean{TypeParameter{Tuple{}}},Euclidean{Tuple{}}}, - x::AbstractVector{<:Number}; + x::AbstractVector{<:Number}, + ::EfficientEstimator; kwargs..., ) return median(x) @@ -551,7 +599,8 @@ end function Statistics.median( ::Union{Euclidean{TypeParameter{Tuple{}}},Euclidean{Tuple{}}}, x::AbstractVector{<:Number}, - w::AbstractWeights; + w::AbstractWeights, + ::EfficientEstimator; kwargs..., ) return median(x, w) diff --git a/src/manifolds/FixedRankMatrices.jl b/src/manifolds/FixedRankMatrices.jl index 97bc66fd18..ddc145ea82 100644 --- a/src/manifolds/FixedRankMatrices.jl +++ b/src/manifolds/FixedRankMatrices.jl @@ -324,15 +324,21 @@ Check whether the tangent [`UMVTVector`](@ref) `X` is from the tangent space of [`FixedRankMatrices`](@ref) `M`, i.e. that `v.U` and `v.Vt` are (columnwise) orthogonal to `x.U` and `x.Vt`, respectively, and its dimensions are consistent with `p` and `X.M`, i.e. correspond to `m`-by-`n` matrices of rank `k`. """ -function check_vector(M::FixedRankMatrices, p::SVDMPoint, X::UMVTVector; kwargs...) +function check_vector( + M::FixedRankMatrices, + p::SVDMPoint, + X::UMVTVector; + atol::Real=sqrt(prod(representation_size(M)) * eps(float(eltype(p.U)))), + kwargs..., +) m, n, k = get_parameter(M.size) - if !isapprox(X.U' * p.U, zeros(k, k); kwargs...) + if !isapprox(X.U' * p.U, zeros(k, k); atol=atol, kwargs...) return DomainError( norm(X.U' * p.U - zeros(k, k)), "The tangent vector $(X) is not a tangent vector to $(p) on $(M) since v.U'x.U is not zero. ", ) end - if !isapprox(X.Vt * p.Vt', zeros(k, k); kwargs...) + if !isapprox(X.Vt * p.Vt', zeros(k, k); atol=atol, kwargs...) return DomainError( norm(X.Vt * p.Vt - zeros(k, k)), "The tangent vector $(X) is not a tangent vector to $(p) on $(M) since v.V'x.V is not zero.", diff --git a/src/manifolds/GeneralUnitaryMatrices.jl b/src/manifolds/GeneralUnitaryMatrices.jl index 5f81197b8d..f6dc9e1df2 100644 --- a/src/manifolds/GeneralUnitaryMatrices.jl +++ b/src/manifolds/GeneralUnitaryMatrices.jl @@ -178,7 +178,7 @@ function cos_angles_4d_rotation_matrix(R) return ((a + b) / 4, (a - b) / 4) end -function default_estimation_method(::GeneralUnitaryMatrices{<:Any,ℝ}, ::typeof(mean)) +function default_approximation_method(::GeneralUnitaryMatrices{<:Any,ℝ}, ::typeof(mean)) return GeodesicInterpolationWithinRadius(π / 2 / √2) end diff --git a/src/manifolds/GeneralizedGrassmann.jl b/src/manifolds/GeneralizedGrassmann.jl index 2fe79dd922..74dd12a183 100644 --- a/src/manifolds/GeneralizedGrassmann.jl +++ b/src/manifolds/GeneralizedGrassmann.jl @@ -280,7 +280,7 @@ Compute the Riemannian [`mean`](@ref mean(M::AbstractManifold, args...)) of `x` """ mean(::GeneralizedGrassmann, ::Any...) -function default_estimation_method(::GeneralizedGrassmann, ::typeof(mean)) +function default_approximation_method(::GeneralizedGrassmann, ::typeof(mean)) return GeodesicInterpolationWithinRadius(π / 4) end diff --git a/src/manifolds/Grassmann.jl b/src/manifolds/Grassmann.jl index c8f23a7612..b3edc40ca8 100644 --- a/src/manifolds/Grassmann.jl +++ b/src/manifolds/Grassmann.jl @@ -171,7 +171,7 @@ Compute the Riemannian [`mean`](@ref mean(M::AbstractManifold, args...)) of `x` """ mean(::Grassmann, ::Any...) -function default_estimation_method(::Grassmann, ::typeof(mean)) +function default_approximation_method(::Grassmann, ::typeof(mean)) return GeodesicInterpolationWithinRadius(π / 4) end diff --git a/src/manifolds/Hyperbolic.jl b/src/manifolds/Hyperbolic.jl index 5fd2e15a10..e6090fea59 100644 --- a/src/manifolds/Hyperbolic.jl +++ b/src/manifolds/Hyperbolic.jl @@ -332,7 +332,7 @@ Compute the Riemannian [`mean`](@ref mean(M::AbstractManifold, args...)) of `x` """ mean(::Hyperbolic, ::Any...) -default_estimation_method(::Hyperbolic, ::typeof(mean)) = CyclicProximalPointEstimation() +default_approximation_method(::Hyperbolic, ::typeof(mean)) = CyclicProximalPointEstimation() @doc raw""" project(M::Hyperbolic, p, X) diff --git a/src/manifolds/HyperbolicHyperboloid.jl b/src/manifolds/HyperbolicHyperboloid.jl index f0b15489b8..b5df6543bb 100644 --- a/src/manifolds/HyperbolicHyperboloid.jl +++ b/src/manifolds/HyperbolicHyperboloid.jl @@ -35,8 +35,14 @@ function check_point(M::Hyperbolic, p; kwargs...) return nothing end -function check_vector(M::Hyperbolic, p, X; kwargs...) - if !isapprox(minkowski_metric(p, X), 0.0; kwargs...) +function check_vector( + M::Hyperbolic, + p, + X::T; + atol::Real=sqrt(prod(representation_size(M))) * eps(real(float(number_eltype(T)))), + kwargs..., +) where {T} + if !isapprox(minkowski_metric(p, X), 0; atol=atol, kwargs...) return DomainError( abs(minkowski_metric(p, X)), "The vector $(X) is not a tangent vector to $(p) on $(M), since it is not orthogonal (with respect to the Minkowski inner product) in the embedding.", diff --git a/src/manifolds/KendallsPreShapeSpace.jl b/src/manifolds/KendallsPreShapeSpace.jl index 3a34719ad9..3b9ef08dc6 100644 --- a/src/manifolds/KendallsPreShapeSpace.jl +++ b/src/manifolds/KendallsPreShapeSpace.jl @@ -37,9 +37,14 @@ representation_size(M::KendallsPreShapeSpace) = get_parameter(M.size) Check whether `p` is a valid point on [`KendallsPreShapeSpace`](@ref), i.e. whether each row has zero mean. Other conditions are checked via embedding in [`ArraySphere`](@ref). """ -function check_point(M::KendallsPreShapeSpace, p; atol=sqrt(eps(eltype(p))), kwargs...) +function check_point( + M::KendallsPreShapeSpace, + p; + atol::Real=sqrt(eps(float(eltype(p)))), + kwargs..., +) for p_row in eachrow(p) - if !isapprox(mean(p_row), 0; atol, kwargs...) + if !isapprox(mean(p_row), 0; atol=atol, kwargs...) return DomainError( mean(p_row), "The point $(p) does not lie on the $(M) since one of the rows does not have zero mean.", @@ -55,9 +60,15 @@ end Check whether `X` is a valid tangent vector on [`KendallsPreShapeSpace`](@ref), i.e. whether each row has zero mean. Other conditions are checked via embedding in [`ArraySphere`](@ref). """ -function check_vector(M::KendallsPreShapeSpace, p, X; atol=sqrt(eps(eltype(X))), kwargs...) +function check_vector( + M::KendallsPreShapeSpace, + p, + X; + atol::Real=sqrt(eps(float(eltype(X)))), + kwargs..., +) for X_row in eachrow(X) - if !isapprox(mean(X_row), 0; atol, kwargs...) + if !isapprox(mean(X_row), 0; atol=atol, kwargs...) return DomainError( mean(X_row), "The vector $(X) is not a tangent vector to $(p) on $(M), since one of the rows does not have zero mean.", diff --git a/src/manifolds/MultinomialDoublyStochastic.jl b/src/manifolds/MultinomialDoublyStochastic.jl index d8dfa1dcac..40a4ef3e03 100644 --- a/src/manifolds/MultinomialDoublyStochastic.jl +++ b/src/manifolds/MultinomialDoublyStochastic.jl @@ -60,10 +60,15 @@ end Checks whether `p` is a valid point on the [`MultinomialDoubleStochastic`](@ref)`(n)` `M`, i.e. is a matrix with positive entries whose rows and columns sum to one. """ -function check_point(M::MultinomialDoubleStochastic, p; kwargs...) +function check_point( + M::MultinomialDoubleStochastic, + p::T; + atol::Real=sqrt(prod(representation_size(M))) * eps(real(float(number_eltype(T)))), + kwargs..., +) where {T} n = get_parameter(M.size)[1] r = sum(p, dims=2) - if !isapprox(norm(r - ones(n, 1)), 0.0; kwargs...) + if !isapprox(norm(r - ones(n, 1)), 0.0; atol=atol, kwargs...) return DomainError( r, "The point $(p) does not lie on $M, since its rows do not sum up to one.", @@ -78,9 +83,15 @@ Checks whether `X` is a valid tangent vector to `p` on the [`MultinomialDoubleSt This means, that `p` is valid, that `X` is of correct dimension and sums to zero along any column or row. """ -function check_vector(M::MultinomialDoubleStochastic, p, X; kwargs...) +function check_vector( + M::MultinomialDoubleStochastic, + p, + X::T; + atol::Real=sqrt(prod(representation_size(M))) * eps(real(float(number_eltype(T)))), + kwargs..., +) where {T} r = sum(X, dims=2) # check for stochastic rows - if !isapprox(norm(r), 0.0; kwargs...) + if !isapprox(norm(r), 0.0; atol=atol, kwargs...) return DomainError( r, "The matrix $(X) is not a tangent vector to $(p) on $(M), since its rows do not sum up to zero.", diff --git a/src/manifolds/ProbabilitySimplex.jl b/src/manifolds/ProbabilitySimplex.jl index 8a00af96c2..96d76fa7b9 100644 --- a/src/manifolds/ProbabilitySimplex.jl +++ b/src/manifolds/ProbabilitySimplex.jl @@ -129,8 +129,14 @@ after [`check_point`](@ref check_point(::ProbabilitySimplex, ::Any))`(M,p)`, `X` has to be of same dimension as `p` and its elements have to sum to one. The tolerance for the last test can be set using the `kwargs...`. """ -function check_vector(M::ProbabilitySimplex, p, X; kwargs...) - if !isapprox(sum(X), 0.0; kwargs...) +function check_vector( + M::ProbabilitySimplex, + p, + X::T; + atol::Real=sqrt(prod(representation_size(M))) * eps(real(float(number_eltype(T)))), + kwargs..., +) where {T} + if !isapprox(sum(X), 0.0; atol=atol, kwargs...) return DomainError( sum(X), "The vector $(X) is not a tangent vector to $(p) on $(M), since its elements do not sum up to 0.", @@ -340,7 +346,7 @@ Compute the Riemannian [`mean`](@ref mean(M::AbstractManifold, args...)) of `x` """ mean(::ProbabilitySimplex, ::Any...) -default_estimation_method(::ProbabilitySimplex, ::typeof(mean)) = GeodesicInterpolation() +default_approximation_method(::ProbabilitySimplex, ::typeof(mean)) = GeodesicInterpolation() function parallel_transport_to!(M::ProbabilitySimplex, Y, p, X, q) n = get_parameter(M.size)[1] diff --git a/src/manifolds/ProjectiveSpace.jl b/src/manifolds/ProjectiveSpace.jl index 392215431b..eeec8e880d 100644 --- a/src/manifolds/ProjectiveSpace.jl +++ b/src/manifolds/ProjectiveSpace.jl @@ -126,8 +126,14 @@ Check whether `X` is a tangent vector in the tangent space of `p` on the tangent space of the embedding and that the Frobenius inner product $⟨p, X⟩_{\mathrm{F}} = 0$. """ -function check_vector(M::AbstractProjectiveSpace, p, X; kwargs...) - if !isapprox(dot(p, X), 0; kwargs...) +function check_vector( + M::AbstractProjectiveSpace, + p, + X::T; + atol::Real=sqrt(prod(representation_size(M))) * eps(real(float(number_eltype(T)))), + kwargs..., +) where {T} + if !isapprox(dot(p, X), 0; atol=atol, kwargs...) return DomainError( dot(p, X), "The vector $(X) is not a tangent vector to $(p) on $(M), since it is not" * @@ -391,7 +397,7 @@ using [`GeodesicInterpolationWithinRadius`](@ref). """ mean(::AbstractProjectiveSpace, ::Any...) -function default_estimation_method(::AbstractProjectiveSpace, ::typeof(mean)) +function default_approximation_method(::AbstractProjectiveSpace, ::typeof(mean)) return GeodesicInterpolationWithinRadius(π / 4) end diff --git a/src/manifolds/Rotations.jl b/src/manifolds/Rotations.jl index b262d0114b..1a1aef5776 100644 --- a/src/manifolds/Rotations.jl +++ b/src/manifolds/Rotations.jl @@ -366,6 +366,12 @@ where ``q=\exp_p d``. The formula simplifies to identity for 2-D rotations. """ parallel_transport_direction(M::Rotations, p, X, d) +function parallel_transport_direction(M::Rotations, p, X, d) + expdhalf = exp(d / 2) + q = exp(M, p, d) + return transpose(q) * p * expdhalf * X * expdhalf +end +parallel_transport_direction(::Rotations{TypeParameter{Tuple{2}}}, p, X, d) = X function parallel_transport_direction!(M::Rotations, Y, p, X, d) expdhalf = exp(d / 2) @@ -375,12 +381,6 @@ end function parallel_transport_direction!(::Rotations{TypeParameter{Tuple{2}}}, Y, p, X, d) return copyto!(Y, X) end -function parallel_transport_direction(M::Rotations, p, X, d) - expdhalf = exp(d / 2) - q = exp(M, p, d) - return transpose(q) * p * expdhalf * X * expdhalf -end -parallel_transport_direction(::Rotations{TypeParameter{Tuple{2}}}, p, X, d) = X function parallel_transport_to!(M::Rotations, Y, p, X, q) d = log(M, p, q) diff --git a/src/manifolds/SPDFixedDeterminant.jl b/src/manifolds/SPDFixedDeterminant.jl index 5efae753c0..d72d3f1f8d 100644 --- a/src/manifolds/SPDFixedDeterminant.jl +++ b/src/manifolds/SPDFixedDeterminant.jl @@ -82,8 +82,14 @@ and additionally fulfill ``\operatorname{tr}(X) = 0``. The tolerance for the trace check of `X` can be set using `kwargs...`, which influences the `isapprox`-check. """ -function check_vector(M::SPDFixedDeterminant, p, X; kwargs...) - if !isapprox(tr(X), 0.0; kwargs...) +function check_vector( + M::SPDFixedDeterminant, + p, + X::T; + atol::Real=sqrt(prod(representation_size(M))) * eps(real(float(number_eltype(T)))), + kwargs..., +) where {T} + if !isapprox(tr(X), 0; atol=atol, kwargs...) return DomainError( tr(X), "The vector $(X) is not a tangent vector to $(p) on $(M), since it does not have a zero trace.", diff --git a/src/manifolds/Spectrahedron.jl b/src/manifolds/Spectrahedron.jl index ed7fb88479..082b6f8fa0 100644 --- a/src/manifolds/Spectrahedron.jl +++ b/src/manifolds/Spectrahedron.jl @@ -85,10 +85,16 @@ and a $X$ has to be a symmetric matrix with trace. The tolerance for the base point check and zero diagonal can be set using the `kwargs...`. Note that symmetry of $X$ holds by construction and is not explicitly checked. """ -function check_vector(M::Spectrahedron, q, Y; kwargs...) +function check_vector( + M::Spectrahedron, + q, + Y::T; + atol::Real=sqrt(prod(representation_size(M))) * eps(real(float(number_eltype(T)))), + kwargs..., +) where {T} X = q * Y' + Y * q' n = tr(X) - if !isapprox(n, 0.0; kwargs...) + if !isapprox(n, 0; atol=atol, kwargs...) return DomainError( n, "The vector $(X) is not a tangent to a point on $(M) (represented py $(q) and $(Y), since its trace is nonzero.", diff --git a/src/manifolds/Sphere.jl b/src/manifolds/Sphere.jl index e2fb758159..73f2cd37c2 100644 --- a/src/manifolds/Sphere.jl +++ b/src/manifolds/Sphere.jl @@ -130,8 +130,14 @@ after [`check_point`](@ref)`(M,p)`, `X` has to be of same dimension as `p` and orthogonal to `p`. The tolerance for the last test can be set using the `kwargs...`. """ -function check_vector(M::AbstractSphere, p, X; kwargs...) - if !isapprox(abs(real(dot(p, X))), 0.0; kwargs...) +function check_vector( + M::AbstractSphere, + p, + X::T; + atol::Real=sqrt(prod(representation_size(M))) * eps(real(float(number_eltype(T)))), + kwargs..., +) where {T} + if !isapprox(abs(real(dot(p, X))), 0; atol=atol, 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.", @@ -413,7 +419,7 @@ Compute the Riemannian [`mean`](@ref mean(M::AbstractManifold, args...)) of `x` """ mean(::AbstractSphere, ::Any...) -function default_estimation_method(::AbstractSphere, ::typeof(mean)) +function default_approximation_method(::AbstractSphere, ::typeof(mean)) return GeodesicInterpolationWithinRadius(π / 2) end diff --git a/src/manifolds/SphereSymmetricMatrices.jl b/src/manifolds/SphereSymmetricMatrices.jl index cdb97c05e5..1c9eecf0f0 100644 --- a/src/manifolds/SphereSymmetricMatrices.jl +++ b/src/manifolds/SphereSymmetricMatrices.jl @@ -35,8 +35,13 @@ i.e. is an `n`-by-`n` symmetric matrix of unit Frobenius norm. The tolerance for the symmetry of `p` can be set using `kwargs...`. """ -function check_point(M::SphereSymmetricMatrices, p; kwargs...) - if !isapprox(norm(p - p'), 0.0; kwargs...) +function check_point( + M::SphereSymmetricMatrices, + p::T; + atol::Real=sqrt(prod(representation_size(M))) * eps(real(float(number_eltype(T)))), + kwargs..., +) where {T} + if !isapprox(norm(p - p'), 0; atol=atol, kwargs...) return DomainError( norm(p - p'), "The point $(p) does not lie on $M, since it is not symmetric.", @@ -54,8 +59,14 @@ of unit Frobenius norm. The tolerance for the symmetry of `p` and `X` can be set using `kwargs...`. """ -function check_vector(M::SphereSymmetricMatrices, p, X; kwargs...) - if !isapprox(norm(X - X'), 0.0; kwargs...) +function check_vector( + M::SphereSymmetricMatrices, + p, + X::T; + atol::Real=sqrt(prod(representation_size(M))) * eps(real(float(number_eltype(T)))), + kwargs..., +) where {T} + if !isapprox(norm(X - X'), 0; atol=atol, kwargs...) return DomainError( norm(X - X'), "The vector $(X) is not a tangent vector to $(p) on $(M), since it is not symmetric.", diff --git a/src/manifolds/Symmetric.jl b/src/manifolds/Symmetric.jl index 3d5a7832bd..989d3dba29 100644 --- a/src/manifolds/Symmetric.jl +++ b/src/manifolds/Symmetric.jl @@ -51,8 +51,8 @@ whether `p` is a symmetric matrix of size `(n,n)` with values from the correspon The tolerance for the symmetry of `p` can be set using `kwargs...`. """ -function check_point(M::SymmetricMatrices{<:Any,𝔽}, p; kwargs...) where {𝔽} - if !isapprox(norm(p - p'), 0.0; kwargs...) +function check_point(M::SymmetricMatrices, p; kwargs...) + if !isapprox(p, p'; kwargs...) return DomainError( norm(p - p'), "The point $(p) does not lie on $M, since it is not symmetric.", @@ -70,8 +70,8 @@ and its values have to be from the correct [`AbstractNumbers`](https://juliamani The tolerance for the symmetry of `X` can be set using `kwargs...`. """ -function check_vector(M::SymmetricMatrices{<:Any,𝔽}, p, X; kwargs...) where {𝔽} - if !isapprox(norm(X - X'), 0.0; kwargs...) +function check_vector(M::SymmetricMatrices, p, X; kwargs...) + if !isapprox(X, X'; kwargs...) return DomainError( norm(X - X'), "The vector $(X) is not a tangent vector to $(p) on $(M), since it is not symmetric.", diff --git a/src/manifolds/SymmetricPositiveDefinite.jl b/src/manifolds/SymmetricPositiveDefinite.jl index 48d367cebd..592f4413eb 100644 --- a/src/manifolds/SymmetricPositiveDefinite.jl +++ b/src/manifolds/SymmetricPositiveDefinite.jl @@ -165,7 +165,7 @@ Lie group. The tolerance for the last test can be set using the `kwargs...`. """ function check_vector(M::SymmetricPositiveDefinite, p, X; kwargs...) - if !isapprox(norm(X - transpose(X)), 0.0; kwargs...) + if !isapprox(X, transpose(X); kwargs...) return DomainError( X, "The vector $(X) is not a tangent to a point on $(M) (represented as an element of the Lie algebra) since its not symmetric.", @@ -289,7 +289,7 @@ Compute the Riemannian [`mean`](@ref mean(M::AbstractManifold, args...)) of `x` """ mean(::SymmetricPositiveDefinite, ::Any) -function default_estimation_method(::SymmetricPositiveDefinite, ::typeof(mean)) +function default_approximation_method(::SymmetricPositiveDefinite, ::typeof(mean)) return GeodesicInterpolation() end diff --git a/src/manifolds/SymmetricPositiveSemidefiniteFixedRank.jl b/src/manifolds/SymmetricPositiveSemidefiniteFixedRank.jl index 0f096829b7..65c774f07c 100644 --- a/src/manifolds/SymmetricPositiveSemidefiniteFixedRank.jl +++ b/src/manifolds/SymmetricPositiveSemidefiniteFixedRank.jl @@ -148,9 +148,15 @@ their distance, if they are not the same, i.e. that $d_{\mathcal M}(p,q) \approx the comparison is performed with the classical `isapprox`. The `kwargs...` are passed on to this accordingly. """ -function _isapprox(M::SymmetricPositiveSemidefiniteFixedRank, p, q; kwargs...) - return isapprox(norm(p - q), 0.0; kwargs...) || - isapprox(distance(M, p, q), 0.0; kwargs...) +function _isapprox( + M::SymmetricPositiveSemidefiniteFixedRank, + p::T, + q; + atol::Real=sqrt(prod(representation_size(M))) * eps(real(float(number_eltype(T)))), + kwargs..., +) where {T} + return isapprox(norm(p - q), 0; atol=atol, kwargs...) || + isapprox(distance(M, p, q), 0; atol=atol, kwargs...) end """ diff --git a/src/manifolds/Symplectic.jl b/src/manifolds/Symplectic.jl index 23cfeca206..67b3492b3f 100644 --- a/src/manifolds/Symplectic.jl +++ b/src/manifolds/Symplectic.jl @@ -212,10 +212,15 @@ Q_{2n} = ```` The tolerance can be set with `kwargs...` (e.g. `atol = 1.0e-14`). """ -function check_point(M::Symplectic, p; kwargs...) +function check_point( + M::Symplectic, + p::T; + atol::Real=sqrt(prod(representation_size(M))) * eps(real(float(number_eltype(T)))), + kwargs..., +) where {T} # Perform check that the matrix lives on the real symplectic manifold: expected_zero = norm(inv(M, p) * p - LinearAlgebra.I) - if !isapprox(expected_zero, zero(eltype(p)); kwargs...) + if !isapprox(expected_zero, 0; atol=atol, kwargs...) return DomainError( expected_zero, ( @@ -245,10 +250,16 @@ The tolerance can be set with `kwargs...` (e.g. `atol = 1.0e-14`). """ check_vector(::Symplectic, ::Any...) -function check_vector(M::Symplectic, p, X; kwargs...) +function check_vector( + M::Symplectic, + p, + X::T; + atol::Real=sqrt(prod(representation_size(M))) * eps(real(float(number_eltype(T)))), + kwargs..., +) where {T} Q = SymplecticMatrix(p, X) tangent_requirement_norm = norm(X' * Q * p + p' * Q * X, 2) - if !isapprox(tangent_requirement_norm, 0.0; kwargs...) + if !isapprox(tangent_requirement_norm, 0; atol=atol, kwargs...) return DomainError( tangent_requirement_norm, ( diff --git a/src/manifolds/SymplecticStiefel.jl b/src/manifolds/SymplecticStiefel.jl index 3a4c1c6ebb..7188e0dafe 100644 --- a/src/manifolds/SymplecticStiefel.jl +++ b/src/manifolds/SymplecticStiefel.jl @@ -103,10 +103,15 @@ Q_{2n} = ```` The tolerance can be set with `kwargs...` (e.g. `atol = 1.0e-14`). """ -function check_point(M::SymplecticStiefel, p; kwargs...) +function check_point( + M::SymplecticStiefel, + p::T; + atol::Real=sqrt(prod(representation_size(M))) * eps(real(float(number_eltype(T)))), + kwargs..., +) where {T} # Perform check that the matrix lives on the real symplectic manifold: expected_zero = norm(inv(M, p) * p - I) - if !isapprox(expected_zero, 0; kwargs...) + if !isapprox(expected_zero, 0; atol=atol, kwargs...) return DomainError( expected_zero, ( @@ -142,14 +147,20 @@ The tolerance can be set with `kwargs...` (e.g. `atol = 1.0e-14`). """ check_vector(::SymplecticStiefel, ::Any...) -function check_vector(M::SymplecticStiefel{<:Any,field}, p, X; kwargs...) where {field} +function check_vector( + M::SymplecticStiefel{S,𝔽}, + p, + X::T; + atol::Real=sqrt(prod(representation_size(M))) * eps(real(float(number_eltype(T)))), + kwargs..., +) where {S,T,𝔽} n, k = get_parameter(M.size) # From Bendokat-Zimmermann: T_pSpSt(2n, 2k) = \{p*H | H^{+} = -H \} H = inv(M, p) * X # ∈ ℝ^{2k × 2k}, should be Hamiltonian. - H_star = inv(Symplectic(2k, field), H) + H_star = inv(Symplectic(2k, 𝔽), H) hamiltonian_identity_norm = norm(H + H_star) - if !isapprox(hamiltonian_identity_norm, 0; kwargs...) + if !isapprox(hamiltonian_identity_norm, 0; atol=atol, kwargs...) return DomainError( hamiltonian_identity_norm, ( diff --git a/src/statistics.jl b/src/statistics.jl index 7a1030e60f..6b6dc83e28 100644 --- a/src/statistics.jl +++ b/src/statistics.jl @@ -1,142 +1,16 @@ """ AbstractEstimationMethod -Abstract type for defining statistical estimation methods. +Deprecated alias for `AbstractApproximationMethod` """ -abstract type AbstractEstimationMethod end - -""" - GradientDescentEstimation <: AbstractEstimationMethod - -Method for estimation using gradient descent. -""" -struct GradientDescentEstimation <: AbstractEstimationMethod end - -""" - CyclicProximalPointEstimation <: AbstractEstimationMethod - -Method for estimation using the cyclic proximal point technique. -""" -struct CyclicProximalPointEstimation <: AbstractEstimationMethod end - -""" - ExtrinsicEstimation <: AbstractEstimationMethod - -Method for estimation in the ambient space and projecting to the manifold. - -For [`mean`](@ref) estimation, [`GeodesicInterpolation`](@ref) is used for mean estimation -in the ambient space. -""" -struct ExtrinsicEstimation <: AbstractEstimationMethod end - -""" - WeiszfeldEstimation <: AbstractEstimationMethod - -Method for estimation using the Weiszfeld algorithm for the [`median`](@ref) -""" -struct WeiszfeldEstimation <: AbstractEstimationMethod end +const AbstractEstimationMethod = AbstractApproximationMethod _unit_weights(n::Int) = StatsBase.UnitWeights{Float64}(n) -@doc raw""" - GeodesicInterpolation <: AbstractEstimationMethod - -Repeated weighted geodesic interpolation method for estimating the Riemannian -center of mass. - -The algorithm proceeds with the following simple online update: - -```math -\begin{aligned} -μ_1 &= x_1\\ -t_k &= \frac{w_k}{\sum_{i=1}^k w_i}\\ -μ_{k} &= γ_{μ_{k-1}}(x_k; t_k), -\end{aligned} -``` - -where $x_k$ are points, $w_k$ are weights, $μ_k$ is the $k$th estimate of the -mean, and $γ_x(y; t)$ is the point at time $t$ along the -[`shortest_geodesic`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/functions.html#ManifoldsBase.shortest_geodesic-Tuple{AbstractManifold,%20Any,%20Any}) -between points $x,y ∈ \mathcal M$. The algorithm -terminates when all $x_k$ have been considered. In the [`Euclidean`](@ref) case, -this exactly computes the weighted mean. - -The algorithm has been shown to converge asymptotically with the sample size for -the following manifolds equipped with their default metrics when all sampled -points are in an open geodesic ball about the mean with corresponding radius -(see [`GeodesicInterpolationWithinRadius`](@ref)): - -* All simply connected complete Riemannian manifolds with non-positive sectional - curvature at radius $∞$ [ChengHoSalehianVemuri:2016](@cite), in particular: - + [`Euclidean`](@ref) - + [`SymmetricPositiveDefinite`](@ref) [HoChengSalehianVemuri:2013](@cite) -* Other manifolds: - + [`Sphere`](@ref): $\frac{π}{2}$ [SalehianEtAl:2015](@cite) - + [`Grassmann`](@ref): $\frac{π}{4}$ [ChakrabortyVemuri:2015](@cite) - + [`Stiefel`](@ref)/[`Rotations`](@ref): $\frac{π}{2 \sqrt 2}$ [ChakrabortyVemuri:2019](@cite) - -For online variance computation, the algorithm additionally uses an analogous -recursion to the weighted Welford algorithm [West:1979](@cite). -""" -struct GeodesicInterpolation <: AbstractEstimationMethod end - -""" - GeodesicInterpolationWithinRadius{T} <: AbstractEstimationMethod - -Estimation of Riemannian center of mass using [`GeodesicInterpolation`](@ref) -with fallback to [`GradientDescentEstimation`](@ref) if any points are outside of a -geodesic ball of specified `radius` around the mean. - -# Constructor - - GeodesicInterpolationWithinRadius(radius) -""" -struct GeodesicInterpolationWithinRadius{T} <: AbstractEstimationMethod - radius::T - - function GeodesicInterpolationWithinRadius(radius::T) where {T} - radius > 0 && return new{T}(radius) - return throw( - DomainError("The radius must be strictly postive, received $(radius)."), - ) - end -end - function Base.show(io::IO, method::GeodesicInterpolationWithinRadius) return print(io, "GeodesicInterpolationWithinRadius($(method.radius))") end -""" - default_estimation_method(M::AbstractManifold, f) - -Specify a default [`AbstractEstimationMethod`](@ref) for an `AbstractManifold` -for a function `f`, e.g. the `median` or the `mean`. - -Note that his function is decorated, so it can inherit from the embedding, for example for the -`IsEmbeddedSubmanifold` trait. -""" -default_estimation_method(M::AbstractManifold, f) - -for mf in [mean, median, cov, var, mean_and_std, mean_and_var] - @eval @trait_function default_estimation_method( - M::AbstractDecoratorManifold, - f::typeof($mf), - ) (no_empty,) - eval( - quote - function default_estimation_method( - ::TraitList{IsEmbeddedSubmanifold}, - M::AbstractDecoratorManifold, - f::typeof($mf), - ) - return default_estimation_method(get_embedding(M), f) - end - end, - ) -end - -@trait_function Statistics.mean(M::AbstractDecoratorManifold, x::AbstractVector) - """ Statistics.cov( M::AbstractManifold, @@ -145,7 +19,7 @@ end tangent_space_covariance_estimator::CovarianceEstimator=SimpleCovariance(; corrected=true, ), - mean_estimation_method::AbstractEstimationMethod=GradientDescentEstimation(), + mean_estimation_method::AbstractApproximationMethod=GradientDescentEstimation(), inverse_retraction_method::AbstractInverseRetractionMethod=default_inverse_retraction_method( M, eltype(x), ), @@ -156,7 +30,7 @@ on a manifold is a rank 2 tensor, the function returns its coefficients in basis the given tangent space basis. See Section 5 of [Pennec:2006](@cite) for details. The mean is calculated using the specified `mean_estimation_method` using -[mean](@ref Statistics.mean(::AbstractManifold, ::AbstractVector, ::AbstractEstimationMethod), +[mean](@ref Statistics.mean(::AbstractManifold, ::AbstractVector, ::AbstractApproximationMethod), and tangent vectors at this mean are calculated using the provided `inverse_retraction_method`. Finally, the covariance matrix in the tangent plane is estimated using the Euclidean space estimator `tangent_space_covariance_estimator`. The type `CovarianceEstimator` is defined @@ -171,7 +45,10 @@ function Statistics.cov( tangent_space_covariance_estimator::CovarianceEstimator=SimpleCovariance(; corrected=true, ), - mean_estimation_method::AbstractEstimationMethod=default_estimation_method(M, cov), + mean_estimation_method::AbstractApproximationMethod=default_approximation_method( + M, + cov, + ), inverse_retraction_method::AbstractInverseRetractionMethod=default_inverse_retraction_method( M, eltype(x), @@ -188,10 +65,16 @@ function Statistics.cov( ) end -function default_estimation_method(::EmptyTrait, ::AbstractDecoratorManifold, ::typeof(cov)) +function default_approximation_method( + ::EmptyTrait, + ::AbstractDecoratorManifold, + ::typeof(cov), +) + return GradientDescentEstimation() +end +function default_approximation_method(::AbstractManifold, ::typeof(cov)) return GradientDescentEstimation() end -default_estimation_method(::AbstractManifold, ::typeof(cov)) = GradientDescentEstimation() @doc raw""" mean(M::AbstractManifold, x::AbstractVector[, w::AbstractWeights]; kwargs...) @@ -209,7 +92,7 @@ In the general case, the [`GradientDescentEstimation`](@ref) is used to compute M::AbstractManifold, x::AbstractVector, [w::AbstractWeights,] - method::AbstractEstimationMethod=default_estimation_method(M); + method::AbstractApproximationMethod=default_approximation_method(M, mean); kwargs..., ) @@ -243,10 +126,25 @@ as the exponential barycenter. The algorithm is further described in[AfsariTronVidal:2013](@cite). """ mean(::AbstractManifold, ::Any...) + +# +# dispatch on method first to allow Euclidean defaults to hit +function Statistics.mean(M::AbstractManifold, x::AbstractVector, kwargs...) + return mean(M, x, default_approximation_method(M, mean, eltype(x)); kwargs...) +end function Statistics.mean( M::AbstractManifold, x::AbstractVector, - method::AbstractEstimationMethod=default_estimation_method(M, mean); + w::AbstractVector, + kwargs..., +) + return mean(M, x, w, default_approximation_method(M, mean, eltype(x)); kwargs...) +end + +function Statistics.mean( + M::AbstractManifold, + x::AbstractVector, + method::AbstractApproximationMethod; kwargs..., ) y = allocate_result(M, mean, x[1]) @@ -256,18 +154,13 @@ function Statistics.mean( M::AbstractManifold, x::AbstractVector, w::AbstractVector, - method::AbstractEstimationMethod=default_estimation_method(M, mean); + method::AbstractApproximationMethod; kwargs..., ) y = allocate_result(M, mean, x[1]) return mean!(M, y, x, w, method; kwargs...) end -function default_estimation_method(::EmptyTrait, ::AbstractManifold, ::typeof(mean)) - return GradientDescentEstimation() -end; -default_estimation_method(::AbstractManifold, ::typeof(mean)) = GradientDescentEstimation(); - @doc raw""" mean!(M::AbstractManifold, y, x::AbstractVector[, w::AbstractWeights]; kwargs...) mean!( @@ -275,7 +168,7 @@ default_estimation_method(::AbstractManifold, ::typeof(mean)) = GradientDescentE y, x::AbstractVector, [w::AbstractWeights,] - method::AbstractEstimationMethod; + method::AbstractApproximationMethod; kwargs..., ) @@ -287,7 +180,7 @@ function Statistics.mean!( M::AbstractManifold, y, x::AbstractVector, - method::AbstractEstimationMethod=default_estimation_method(M, mean); + method::AbstractApproximationMethod=default_approximation_method(M, mean); kwargs..., ) w = _unit_weights(length(x)) @@ -497,8 +390,6 @@ end Estimate the Riemannian center of mass of `x` using [`ExtrinsicEstimation`](@ref), i.e. by computing the mean in the embedding and projecting the result back. -You can specify an `extrinsic_method` to specify which mean estimation method to use in the embedding, -which defaults to [`GeodesicInterpolation`](@ref). See [`mean`](@ref mean(::AbstractManifold, ::AbstractVector, ::AbstractVector, ::GeodesicInterpolation)) for a description of the remaining `kwargs`. @@ -515,26 +406,37 @@ function Statistics.mean!( y, x::AbstractVector, w::AbstractVector, - ::ExtrinsicEstimation; - extrinsic_method::AbstractEstimationMethod=default_estimation_method( - get_embedding(M), - mean, - ), + e::ExtrinsicEstimation; + extrinsic_method::Union{AbstractEstimationMethod,Nothing}=nothing, kwargs..., ) + if !isnothing(extrinsic_method) + Base.depwarn( + "The Keyword Argument `extrinsic_method` is deprecated use `ExtrinsicEstimators` field instead", + :mean!, + ) + e = ExtrinsicEstimation(extrinsic_method) + end embedded_x = map(p -> embed(M, p), x) - embedded_y = mean(get_embedding(M), embedded_x, w, extrinsic_method; kwargs...) + embedded_y = mean(get_embedding(M), embedded_x, w, e.extrinsic_estimation; kwargs...) project!(M, y, embedded_y) return y end +function default_approximation_method(::EmptyTrait, ::AbstractManifold, ::typeof(mean)) + return GradientDescentEstimation() +end; +function default_approximation_method(::AbstractManifold, ::typeof(mean)) + return GradientDescentEstimation() +end; + @doc raw""" median(M::AbstractManifold, x::AbstractVector[, w::AbstractWeights]; kwargs...) median( M::AbstractManifold, x::AbstractVector, [w::AbstractWeights,] - method::AbstractEstimationMethod; + method::AbstractApproximationMethod; kwargs..., ) @@ -553,14 +455,14 @@ Compute the median using the specified `method`. """ Statistics.median(::AbstractManifold, ::Any...) -function default_estimation_method( +function default_approximation_method( ::EmptyTrait, ::AbstractDecoratorManifold, ::typeof(median), ) return CyclicProximalPointEstimation() end -function default_estimation_method(::AbstractManifold, ::typeof(median)) +function default_approximation_method(::AbstractManifold, ::typeof(median)) return CyclicProximalPointEstimation() end @@ -605,14 +507,11 @@ Statistics.median( x::AbstractVector, [w::AbstractWeights,] method::ExtrinsicEstimation; - extrinsic_method = CyclicProximalPointEstimation(), kwargs..., ) Estimate the median of `x` using [`ExtrinsicEstimation`](@ref), i.e. by computing the median in the embedding and projecting the result back. -You can specify an `extrinsic_method` to specify which median estimation method to use in -the embedding, which defaults to [`CyclicProximalPointEstimation`](@ref). See [`median`](@ref median(::AbstractManifold, ::AbstractVector, ::AbstractVector, ::CyclicProximalPointEstimation)) for a description of `kwargs`. @@ -680,10 +579,24 @@ Statistics.median( ::WeiszfeldEstimation, ) +# +# dispatch on the method first before allocating to allow Euclidean defaults to hit +function Statistics.median(M::AbstractManifold, x::AbstractVector; kwargs...) + return median(M, x, default_approximation_method(M, median, eltype(x))) +end function Statistics.median( M::AbstractManifold, x::AbstractVector, - method::AbstractEstimationMethod=default_estimation_method(M, median); + w::AbstractVector; + kwargs..., +) + return median(M, x, w, default_approximation_method(M, median, eltype(x))) +end + +function Statistics.median( + M::AbstractManifold, + x::AbstractVector, + method::AbstractApproximationMethod; kwargs..., ) y = allocate_result(M, median, x[1]) @@ -693,7 +606,7 @@ function Statistics.median( M::AbstractManifold, x::AbstractVector, w::AbstractVector, - method::AbstractEstimationMethod=default_estimation_method(M, median); + method::AbstractApproximationMethod; kwargs..., ) y = allocate_result(M, median, x[1]) @@ -707,7 +620,7 @@ end y, x::AbstractVector, [w::AbstractWeights,] - method::AbstractEstimationMethod; + method::AbstractApproximationMethod; kwargs..., ) @@ -718,7 +631,7 @@ function Statistics.median!( M::AbstractManifold, q, x::AbstractVector, - method::AbstractEstimationMethod=default_estimation_method(M, median); + method::AbstractApproximationMethod=default_approximation_method(M, median, eltype(x)); kwargs..., ) w = _unit_weights(length(x)) @@ -771,15 +684,19 @@ function Statistics.median!( y, x::AbstractVector, w::AbstractVector, - ::ExtrinsicEstimation; - extrinsic_method::AbstractEstimationMethod=default_estimation_method( - get_embedding(M), - median, - ), + e::ExtrinsicEstimation; + extrinsic_method=nothing, kwargs..., ) + if !isnothing(extrinsic_method) + Base.depwarn( + "The Keyword Argument `extrinsic_method` is deprecated use `ExtrinsicEstimators` field instead", + :median!, + ) + e = ExtrinsicEstimation(extrinsic_method) + end embedded_x = map(p -> embed(M, p), x) - embedded_y = median(get_embedding(M), embedded_x, w, extrinsic_method; kwargs...) + embedded_y = median(get_embedding(M), embedded_x, w, e.extrinsic_estimation; kwargs...) project!(M, y, embedded_y) return y end @@ -905,7 +822,7 @@ simultaneously. See those functions for a description of the arguments. M::AbstractManifold, x::AbstractVector [w::AbstractWeights,] - method::AbstractEstimationMethod; + method::AbstractApproximationMethod; kwargs..., ) -> (mean, var) @@ -918,7 +835,7 @@ function StatsBase.mean_and_var( M::AbstractManifold, x::AbstractVector, w::AbstractWeights, - method::AbstractEstimationMethod=default_estimation_method(M, mean); + method::AbstractApproximationMethod=default_approximation_method(M, mean, eltype(x)); corrected=false, kwargs..., ) @@ -929,7 +846,11 @@ end function StatsBase.mean_and_var( M::AbstractManifold, x::AbstractVector, - method::AbstractEstimationMethod=default_estimation_method(M, mean_and_var); + method::AbstractApproximationMethod=default_approximation_method( + M, + mean_and_var, + eltype(x), + ); corrected=true, kwargs..., ) @@ -937,15 +858,15 @@ function StatsBase.mean_and_var( w = _unit_weights(n) return mean_and_var(M, x, w, method; corrected=corrected, kwargs...) end -function default_estimation_method( +function default_approximation_method( ::EmptyTrait, M::AbstractDecoratorManifold, ::typeof(mean_and_var), ) - return default_estimation_method(M, mean) + return default_approximation_method(M, mean) end -function default_estimation_method(M::AbstractManifold, ::typeof(mean_and_var)) - return default_estimation_method(M, mean) +function default_approximation_method(M::AbstractManifold, ::typeof(mean_and_var)) + return default_approximation_method(M, mean) end @doc raw""" @@ -1077,7 +998,7 @@ Compute the [`mean`](@ref mean(::AbstractManifold, args...)) and the standard de M::AbstractManifold, x::AbstractVector [w::AbstractWeights,] - method::AbstractEstimationMethod; + method::AbstractApproximationMethod; kwargs..., ) -> (mean, var) @@ -1089,8 +1010,8 @@ function StatsBase.mean_and_std(M::AbstractManifold, args...; kwargs...) m, v = mean_and_var(M, args...; kwargs...) return m, sqrt(v) end -function default_estimation_method(M::AbstractManifold, ::typeof(mean_and_std)) - return default_estimation_method(M, mean) +function default_approximation_method(M::AbstractManifold, ::typeof(mean_and_std)) + return default_approximation_method(M, mean) end """ @@ -1154,3 +1075,24 @@ function StatsBase.kurtosis(M::AbstractManifold, x::AbstractVector, args...) w = _unit_weights(length(x)) return kurtosis(M, x, w, args...) end + +# +# decorate default method for a few functions +for mf in [mean, median, cov, var, mean_and_std, mean_and_var] + @eval @trait_function default_approximation_method( + M::AbstractDecoratorManifold, + f::typeof($mf), + ) (no_empty,) + eval( + quote + function default_approximation_method( + ::TraitList{IsEmbeddedSubmanifold}, + M::AbstractDecoratorManifold, + f::typeof($mf), + ) + return default_approximation_method(get_embedding(M), f) + end + end, + ) +end +@trait_function Statistics.mean(M::AbstractDecoratorManifold, x::AbstractVector) diff --git a/test/manifolds/circle.jl b/test/manifolds/circle.jl index 07d3eddb6d..46bfe9d226 100644 --- a/test/manifolds/circle.jl +++ b/test/manifolds/circle.jl @@ -288,4 +288,10 @@ using Manifolds: TFVector, CoTFVector ) end end + @testset "Mixed array dimensions for exp" begin + M = Circle() + p = fill(0.0) + exp!(M, p, p, [1.0], 2.0) + @test p ≈ fill(2.0) + end end diff --git a/test/manifolds/embedded_torus.jl b/test/manifolds/embedded_torus.jl index a5ccf63e36..b8720d10d5 100644 --- a/test/manifolds/embedded_torus.jl +++ b/test/manifolds/embedded_torus.jl @@ -91,7 +91,7 @@ using BoundaryValueDiffEq a2 = [-0.5, 0.3] sol_log = Manifolds.solve_chart_log_bvp(M, p0x, a2, A, (0, 0)) @test sol_log(0.0)[1:2] ≈ p0x - @test sol_log(1.0)[1:2] ≈ a2 + @test sol_log(1.0)[1:2] ≈ a2 atol = 1e-7 # a test randomly failed here on Julia 1.6 once for no clear reason? # so I bumped tolerance considerably bvp_atol = VERSION < v"1.7" ? 2e-3 : 1e-15 diff --git a/test/manifolds/rotations.jl b/test/manifolds/rotations.jl index 3fa56e491b..4a276b8def 100644 --- a/test/manifolds/rotations.jl +++ b/test/manifolds/rotations.jl @@ -286,4 +286,33 @@ include("../utils.jl") @test X isa Matrix{Float64} @test X == fill(0.0, 1, 1) end + + @testset "Specializations" begin + M = Rotations(2) + p = Matrix{Float64}(I, 2, 2) + X = [0.0 3.0; -3.0 0.0] + @test parallel_transport_direction(M, p, X, X) === X + + M = Rotations(3) + p = @SMatrix [ + -0.5908399013383766 -0.6241917041179139 0.5111681988316876 + -0.7261666986267721 0.13535732881097293 -0.6740625485388226 + 0.35155388888753836 -0.7694563730631729 -0.5332417398896261 + ] + X = @SMatrix [ + 0.0 -0.30777760628130063 0.5499897386953444 + 0.30777760628130063 0.0 -0.32059980100053004 + -0.5499897386953444 0.32059980100053004 0.0 + ] + d = @SMatrix [ + 0.0 -0.4821890003925358 -0.3513148535122392 + 0.4821890003925358 0.0 0.37956770358148356 + 0.3513148535122392 -0.37956770358148356 0.0 + ] + @test parallel_transport_direction(M, p, X, d) ≈ [ + 0.0 -0.3258778314599828 0.3903114578816008 + 0.32587783145998306 0.0 -0.49138641089195584 + -0.3903114578816011 0.4913864108919558 0.0 + ] + end end diff --git a/test/statistics.jl b/test/statistics.jl index c4356d3788..160537bc00 100644 --- a/test/statistics.jl +++ b/test/statistics.jl @@ -12,13 +12,14 @@ import ManifoldsBase: base_manifold, get_embedding using Manifolds: - AbstractEstimationMethod, + AbstractApproximationMethod, CyclicProximalPointEstimation, GeodesicInterpolation, GeodesicInterpolationWithinRadius, GradientDescentEstimation, WeiszfeldEstimation -import Manifolds: mean, mean!, median, median!, var, mean_and_var, default_estimation_method +import Manifolds: + mean, mean!, median, median!, var, mean_and_var, default_approximation_method struct TestStatsSphere{N} <: AbstractManifold{ℝ} end TestStatsSphere(N) = TestStatsSphere{N}() @@ -114,7 +115,7 @@ function test_mean(M, x, yexp=nothing, method...; kwargs...) y, x, pweights(ones(n + 1)), - Manifolds.default_estimation_method(M, mean); + Manifolds.default_approximation_method(M, mean); kwargs..., ) end @@ -125,14 +126,14 @@ function test_median( M, x, yexp=nothing; - method::Union{Nothing,AbstractEstimationMethod}=nothing, + method::Union{Nothing,AbstractApproximationMethod}=nothing, kwargs..., ) @testset "median unweighted$(!isnothing(method) ? " ($method)" : "")" begin y = isnothing(method) ? median(M, x; kwargs...) : median(M, x, method; kwargs...) @test is_point(M, y; atol=10^-9) if yexp !== nothing - @test isapprox(M, y, yexp; atol=10^-5) + @test isapprox(M, y, yexp; atol=5 * 10^-5) end end @@ -301,7 +302,7 @@ function test_moments(M, x) end struct TestStatsOverload1 <: AbstractManifold{ℝ} end -struct TestStatsMethod1 <: AbstractEstimationMethod end +struct TestStatsMethod1 <: AbstractApproximationMethod end function mean!( ::TestStatsOverload1, @@ -401,8 +402,8 @@ end @test std(M, x, w) == 2.0 @test std(M, x, w, 2) == 2.0 - @test Manifolds.default_estimation_method(M, mean_and_std) == - Manifolds.default_estimation_method(M, mean) + @test Manifolds.default_approximation_method(M, mean_and_std) == + Manifolds.default_approximation_method(M, mean) @test mean_and_var(M, x, TestStatsMethod1()) == ([5.0], 16) @test mean_and_var(M, x, w, TestStatsMethod1()) == ([5.0], 9) @test mean_and_std(M, x, TestStatsMethod1()) == ([5.0], 4.0) @@ -547,6 +548,9 @@ end test_var(M, x) test_std(M, x) test_moments(M, x) + y = copy(x[1]) + mean!(M, y, x) + @test y == mean(x) end end end @@ -786,7 +790,7 @@ end x = [normalize(randn(rng, 3)) for _ in 1:10] w = pweights([rand(rng) for _ in 1:length(x)]) m = normalize(mean(reduce(hcat, x), w; dims=2)[:, 1]) - mg = mean(S, x, w, ExtrinsicEstimation()) + mg = mean(S, x, w, ExtrinsicEstimation(EfficientEstimator())) @test isapprox(S, m, mg) end @@ -796,12 +800,12 @@ end x = [normalize(randn(rng, 3)) for _ in 1:10] w = pweights([rand(rng) for _ in 1:length(x)]) m = normalize(median(Euclidean(3), x, w)) - mg = median(S, x, w, ExtrinsicEstimation()) + mg = median(S, x, w, ExtrinsicEstimation(CyclicProximalPointEstimation())) @test isapprox(S, m, mg) end @testset "Covariance Default" begin - @test default_estimation_method(TestStatsSphere{2}(), cov) == + @test default_approximation_method(TestStatsSphere{2}(), cov) == GradientDescentEstimation() end diff --git a/test/test_deprecated.jl b/test/test_deprecated.jl index 5791392758..1fb4c592f5 100644 --- a/test/test_deprecated.jl +++ b/test/test_deprecated.jl @@ -1,3 +1,33 @@ -using Manifolds, ManifoldsBase, Test +using Manifolds, ManifoldsBase, Random, Test -@testset "Deprecation tests" begin end +using StatsBase: AbstractWeights, pweights +using Random: GLOBAL_RNG, seed! + +@testset "Deprecation tests" begin + @testset "Depreacte extrinsic_method= keyword" begin + rng = MersenneTwister(47) + S = Sphere(2) + x = [normalize(randn(rng, 3)) for _ in 1:10] + w = pweights([rand(rng) for _ in 1:length(x)]) + mg1 = mean(S, x, w, ExtrinsicEstimation(EfficientEstimator())) + # Statistics 414-418, depcreatce former extrinsic_method keyword + mg2 = mean( + S, + x, + w, + ExtrinsicEstimation(EfficientEstimator()); + extrinsic_method=EfficientEstimator(), + ) + @test isapprox(S, mg1, mg2) + mg3 = median(S, x, w, ExtrinsicEstimation(CyclicProximalPointEstimation())) + # Statistics 692-696, depcreatce former extrinsic_method keyword + mg4 = median( + S, + x, + w, + ExtrinsicEstimation(CyclicProximalPointEstimation()); + extrinsic_method=CyclicProximalPointEstimation(), + ) + @test isapprox(S, mg3, mg4) + end +end From 631e8b45669c854e6275a10c64a4a7b53ad5c9fe Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Mon, 25 Dec 2023 18:07:32 +0100 Subject: [PATCH 2/5] Fix documentation on master (#694) * Fix most of the docs. * Update news entry. --- NEWS.md | 2 +- docs/Project.toml | 2 +- docs/src/features/statistics.md | 2 +- docs/src/manifolds/fixedrankmatrices.md | 2 +- docs/src/manifolds/grassmann.md | 4 ++-- docs/src/manifolds/hyperbolic.md | 2 +- docs/src/manifolds/probabilitysimplex.md | 2 +- docs/src/manifolds/rotations.md | 2 +- docs/src/manifolds/sphere.md | 2 +- docs/src/manifolds/stiefel.md | 2 +- docs/src/manifolds/symplectic.md | 2 +- docs/src/manifolds/symplecticstiefel.md | 2 +- docs/src/references.bib | 8 ++++---- src/manifolds/GeneralUnitaryMatrices.jl | 2 +- src/manifolds/GeneralizedGrassmann.jl | 2 +- src/manifolds/Grassmann.jl | 2 +- src/manifolds/Hyperbolic.jl | 2 +- src/manifolds/ProbabilitySimplex.jl | 2 +- src/manifolds/ProjectiveSpace.jl | 2 +- src/manifolds/Sphere.jl | 2 +- src/manifolds/SymmetricPositiveDefinite.jl | 2 +- src/statistics.jl | 22 +++++++++++----------- tutorials/getstarted.qmd | 2 +- tutorials/integration.qmd | 2 +- 24 files changed, 38 insertions(+), 38 deletions(-) diff --git a/NEWS.md b/NEWS.md index 89c862fe12..123ea7737b 100644 --- a/NEWS.md +++ b/NEWS.md @@ -5,7 +5,7 @@ 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.9] – 2023-12-24 +## [0.9.9] – 2023-12-25 ### Fixed diff --git a/docs/Project.toml b/docs/Project.toml index b9d1637250..b7ce0d6e09 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -25,7 +25,7 @@ CondaPkg = "0.2" DiffEqCallbacks = "2" Distributions = "0.22.6, 0.23, 0.24, 0.25" Documenter = "1" -DocumenterCitations = "1.2.1" +DocumenterCitations = "1.3" FiniteDifferences = "0.12" Graphs = "1.4" HybridArrays = "0.4" diff --git a/docs/src/features/statistics.md b/docs/src/features/statistics.md index 51df528995..716b5172e4 100644 --- a/docs/src/features/statistics.md +++ b/docs/src/features/statistics.md @@ -9,6 +9,6 @@ Order = [:type, :function] ## Literature ```@bibliography -Pages = ["features/statistics.md"] +Pages = ["statistics.md"] Canonical=false ``` \ No newline at end of file diff --git a/docs/src/manifolds/fixedrankmatrices.md b/docs/src/manifolds/fixedrankmatrices.md index c1fb92ab7a..e2ff785e00 100644 --- a/docs/src/manifolds/fixedrankmatrices.md +++ b/docs/src/manifolds/fixedrankmatrices.md @@ -9,6 +9,6 @@ Order = [:type, :function] ## Literature ```@bibliography -Pages = ["manifolds/fixedrankmatrices.md"] +Pages = ["fixedrankmatrices.md"] Canonical=false ``` \ No newline at end of file diff --git a/docs/src/manifolds/grassmann.md b/docs/src/manifolds/grassmann.md index d0e964592f..1fb337c887 100644 --- a/docs/src/manifolds/grassmann.md +++ b/docs/src/manifolds/grassmann.md @@ -10,7 +10,7 @@ Order = [:type,:function] ```@autodocs Modules = [Manifolds] -Pages = ["manifolds/GrassmannStiefel.jl"] +Pages = ["GrassmannStiefel.jl"] Order = [:type,:function] ``` @@ -26,6 +26,6 @@ Order = [:type,:function] ## Literature ```@bibliography -Pages = ["manifolds/grassmann.md"] +Pages = ["grassmann.md"] Canonical=false ``` \ No newline at end of file diff --git a/docs/src/manifolds/hyperbolic.md b/docs/src/manifolds/hyperbolic.md index e3cb4bc089..c9976bbf4f 100644 --- a/docs/src/manifolds/hyperbolic.md +++ b/docs/src/manifolds/hyperbolic.md @@ -180,6 +180,6 @@ plot!(scene, M, [pts[2], origin], geodesic_interpolation=100) ## Literature ```@bibliography -Pages = ["manifolds/hyperbolic.md"] +Pages = ["hyperbolic.md"] Canonical=false ``` \ No newline at end of file diff --git a/docs/src/manifolds/probabilitysimplex.md b/docs/src/manifolds/probabilitysimplex.md index 182b838393..bc7b76c71a 100644 --- a/docs/src/manifolds/probabilitysimplex.md +++ b/docs/src/manifolds/probabilitysimplex.md @@ -43,6 +43,6 @@ Public=false ## Literature ```@bibliography -Pages = ["manifolds/probabilitysimplex.md"] +Pages = ["probabilitysimplex.md"] Canonical=false ``` \ No newline at end of file diff --git a/docs/src/manifolds/rotations.md b/docs/src/manifolds/rotations.md index 5eafa1f249..0458a7e58c 100644 --- a/docs/src/manifolds/rotations.md +++ b/docs/src/manifolds/rotations.md @@ -39,6 +39,6 @@ Order = [:type, :function] ## Literature ```@bibliography -Pages = ["manifolds/rotations.md"] +Pages = ["rotations.md"] Canonical=false ``` diff --git a/docs/src/manifolds/sphere.md b/docs/src/manifolds/sphere.md index 3153ffddc5..c4ae99589e 100644 --- a/docs/src/manifolds/sphere.md +++ b/docs/src/manifolds/sphere.md @@ -66,6 +66,6 @@ plot!(scene, M, pts2, vecs; wireframe = false, linewidth=1.5) ## Literature ```@bibliography -Pages = ["manifolds/sphere.md"] +Pages = ["sphere.md"] Canonical=false ``` \ No newline at end of file diff --git a/docs/src/manifolds/stiefel.md b/docs/src/manifolds/stiefel.md index ca988418e2..66019c36c1 100644 --- a/docs/src/manifolds/stiefel.md +++ b/docs/src/manifolds/stiefel.md @@ -67,6 +67,6 @@ Private = true ## Literature ```@bibliography -Pages = ["manifolds/stiefel.md"] +Pages = ["stiefel.md"] Canonical=false ``` \ No newline at end of file diff --git a/docs/src/manifolds/symplectic.md b/docs/src/manifolds/symplectic.md index 7f435c10b3..3673e6be2b 100644 --- a/docs/src/manifolds/symplectic.md +++ b/docs/src/manifolds/symplectic.md @@ -40,6 +40,6 @@ Order = [:type, :function] ## Literature ```@bibliography -Pages = ["manifolds/symplectic.md"] +Pages = ["symplectic.md"] Canonical=false ``` \ No newline at end of file diff --git a/docs/src/manifolds/symplecticstiefel.md b/docs/src/manifolds/symplecticstiefel.md index 2fda3518e8..502f97e211 100644 --- a/docs/src/manifolds/symplecticstiefel.md +++ b/docs/src/manifolds/symplecticstiefel.md @@ -31,6 +31,6 @@ Order = [:type, :function] ## Literature ```@bibliography -Pages = ["manifolds/symplecticstiefel.md"] +Pages = ["symplecticstiefel.md"] Canonical=false ``` \ No newline at end of file diff --git a/docs/src/references.bib b/docs/src/references.bib index 0b43a860f6..98a24f12f6 100644 --- a/docs/src/references.bib +++ b/docs/src/references.bib @@ -10,7 +10,7 @@ @book{AbsilMahonySepulchre:2008, AUTHOR = {Absil, P.-A. and Mahony, R. and Sepulchre, R.}, DOI = {10.1515/9781400830244}, - NOTE = {available online at [press.princeton.edu/chapters/absil/](http://press.princeton.edu/chapters/absil/)}, + NOTE = {available online at \href{http://press.princeton.edu/chapters/absil/}{press.princeton.edu/chapters/absil/}}, PUBLISHER = {Princeton University Press}, TITLE = {Optimization Algorithms on Matrix Manifolds}, YEAR = {2008}, @@ -108,7 +108,6 @@ @article{Bacak:2014 EPRINT = {1210.2145}, EPRINTTYPE = {arXiv}, JOURNAL = {SIAM Journal on Optimization}, - NOTE = {arXiv: [1210.2145](https://arxiv.org/abs/1210.2145)}, NUMBER = {3}, PAGES = {1542--1566}, TITLE = {Computing medians and means in Hadamard spaces}, @@ -117,9 +116,10 @@ @article{Bacak:2014 } @article{BendokatZimmermann:2021, AUTHOR = {Bendokat, Thomas and Zimmermann, Ralf}, + EPRINT = {2108.12447}, + EPRINTTYPE = {arXiv}, Journal = {arXiv Preprint, 2108.12447}, TITLE = {The real symplectic Stiefel and Grassmann manifolds: metrics, geodesics and applications}, - URL = {https://arxiv.org/abs/2108.12447}, YEAR = {2021} } @article{BendokatZimmermannAbsil:2020, @@ -588,7 +588,7 @@ @article{MassartAbsil:2020 DOI = {10.1137/18m1231389}, YEAR = {2020}, VOLUME = {41}, - NOTE = {Preprint: [sites.uclouvain.be/absil/2018.06](https://sites.uclouvain.be/absil/2018.06)}, + NOTE = {Preprint: \href{https://sites.uclouvain.be/absil/2018.06}{sites.uclouvain.be/absil/2018.06}}, NUMBER = {1}, PAGES = {171--198}, AUTHOR = {Estelle Massart and P.-A. Absil}, diff --git a/src/manifolds/GeneralUnitaryMatrices.jl b/src/manifolds/GeneralUnitaryMatrices.jl index f6dc9e1df2..6ff6f9068c 100644 --- a/src/manifolds/GeneralUnitaryMatrices.jl +++ b/src/manifolds/GeneralUnitaryMatrices.jl @@ -946,7 +946,7 @@ end ) Compute the Riemannian [`mean`](@ref mean(M::AbstractManifold, args...)) of `x` using -[`GeodesicInterpolationWithinRadius`](@ref). +[`GeodesicInterpolationWithinRadius`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/functions/#ManifoldsBase.GeodesicInterpolationWithinRadius). """ mean(::GeneralUnitaryMatrices{<:Any,ℝ}, ::Any) diff --git a/src/manifolds/GeneralizedGrassmann.jl b/src/manifolds/GeneralizedGrassmann.jl index 74dd12a183..7fd9c26dff 100644 --- a/src/manifolds/GeneralizedGrassmann.jl +++ b/src/manifolds/GeneralizedGrassmann.jl @@ -276,7 +276,7 @@ end ) Compute the Riemannian [`mean`](@ref mean(M::AbstractManifold, args...)) of `x` using -[`GeodesicInterpolationWithinRadius`](@ref). +[`GeodesicInterpolationWithinRadius`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/functions/#ManifoldsBase.GeodesicInterpolationWithinRadius). """ mean(::GeneralizedGrassmann, ::Any...) diff --git a/src/manifolds/Grassmann.jl b/src/manifolds/Grassmann.jl index b3edc40ca8..a7adb80612 100644 --- a/src/manifolds/Grassmann.jl +++ b/src/manifolds/Grassmann.jl @@ -167,7 +167,7 @@ end ) Compute the Riemannian [`mean`](@ref mean(M::AbstractManifold, args...)) of `x` using -[`GeodesicInterpolationWithinRadius`](@ref). +[`GeodesicInterpolationWithinRadius`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/functions/#ManifoldsBase.GeodesicInterpolationWithinRadius). """ mean(::Grassmann, ::Any...) diff --git a/src/manifolds/Hyperbolic.jl b/src/manifolds/Hyperbolic.jl index e6090fea59..3d40bbcd4c 100644 --- a/src/manifolds/Hyperbolic.jl +++ b/src/manifolds/Hyperbolic.jl @@ -328,7 +328,7 @@ manifold_volume(::Hyperbolic) = Inf ) Compute the Riemannian [`mean`](@ref mean(M::AbstractManifold, args...)) of `x` on the -[`Hyperbolic`](@ref) space using [`CyclicProximalPointEstimation`](@ref). +[`Hyperbolic`](@ref) space using [`CyclicProximalPointEstimation`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/functions/#ManifoldsBase.CyclicProximalPointEstimation). """ mean(::Hyperbolic, ::Any...) diff --git a/src/manifolds/ProbabilitySimplex.jl b/src/manifolds/ProbabilitySimplex.jl index 96d76fa7b9..7a7983314d 100644 --- a/src/manifolds/ProbabilitySimplex.jl +++ b/src/manifolds/ProbabilitySimplex.jl @@ -342,7 +342,7 @@ end ) Compute the Riemannian [`mean`](@ref mean(M::AbstractManifold, args...)) of `x` using -[`GeodesicInterpolation`](@ref). +[`GeodesicInterpolation`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/functions/#ManifoldsBase.GeodesicInterpolation). """ mean(::ProbabilitySimplex, ::Any...) diff --git a/src/manifolds/ProjectiveSpace.jl b/src/manifolds/ProjectiveSpace.jl index eeec8e880d..7135eb4df2 100644 --- a/src/manifolds/ProjectiveSpace.jl +++ b/src/manifolds/ProjectiveSpace.jl @@ -393,7 +393,7 @@ end ) Compute the Riemannian [`mean`](@ref mean(M::AbstractManifold, args...)) of points in vector `x` -using [`GeodesicInterpolationWithinRadius`](@ref). +using [`GeodesicInterpolationWithinRadius`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/functions/#ManifoldsBase.GeodesicInterpolationWithinRadius). """ mean(::AbstractProjectiveSpace, ::Any...) diff --git a/src/manifolds/Sphere.jl b/src/manifolds/Sphere.jl index 73f2cd37c2..b39be325f1 100644 --- a/src/manifolds/Sphere.jl +++ b/src/manifolds/Sphere.jl @@ -415,7 +415,7 @@ end ) Compute the Riemannian [`mean`](@ref mean(M::AbstractManifold, args...)) of `x` using -[`GeodesicInterpolationWithinRadius`](@ref). +[`GeodesicInterpolationWithinRadius`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/functions/#ManifoldsBase.GeodesicInterpolationWithinRadius). """ mean(::AbstractSphere, ::Any...) diff --git a/src/manifolds/SymmetricPositiveDefinite.jl b/src/manifolds/SymmetricPositiveDefinite.jl index 592f4413eb..13b58b171b 100644 --- a/src/manifolds/SymmetricPositiveDefinite.jl +++ b/src/manifolds/SymmetricPositiveDefinite.jl @@ -285,7 +285,7 @@ end ) Compute the Riemannian [`mean`](@ref mean(M::AbstractManifold, args...)) of `x` using -[`GeodesicInterpolation`](@ref). +[`GeodesicInterpolation`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/functions/#ManifoldsBase.GeodesicInterpolation). """ mean(::SymmetricPositiveDefinite, ::Any) diff --git a/src/statistics.jl b/src/statistics.jl index 6b6dc83e28..e4cad6c6fd 100644 --- a/src/statistics.jl +++ b/src/statistics.jl @@ -87,7 +87,7 @@ as the point that satisfies the minimizer ```` where $\mathrm{d}_{\mathcal M}$ denotes the Riemannian [`distance`](@ref). -In the general case, the [`GradientDescentEstimation`](@ref) is used to compute the mean. +In the general case, the [`GradientDescentEstimation`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/functions/#ManifoldsBase.GradientDescentEstimation) is used to compute the mean. mean( M::AbstractManifold, x::AbstractVector, @@ -110,7 +110,7 @@ Compute the mean using the specified `method`. kwargs..., ) -Compute the mean using the gradient descent scheme [`GradientDescentEstimation`](@ref). +Compute the mean using the gradient descent scheme [`GradientDescentEstimation`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/functions/#ManifoldsBase.GradientDescentEstimation). Optionally, provide `p0`, the starting point (by default set to the first data point). `stop_iter` denotes the maximal number of iterations to perform and the @@ -243,7 +243,7 @@ end Estimate the Riemannian center of mass of `x` in an online fashion using repeated weighted geodesic interpolation. See -[`GeodesicInterpolation`](@ref) for details. +[`GeodesicInterpolation`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/functions/#ManifoldsBase.GeodesicInterpolation) for details. If `shuffle_rng` is provided, it is used to shuffle the order in which the points are considered for computing the mean. @@ -305,7 +305,7 @@ end ) Estimate the Riemannian center of mass of `x` using -[`GeodesicInterpolationWithinRadius`](@ref). +[`GeodesicInterpolationWithinRadius`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/functions/#ManifoldsBase.GeodesicInterpolationWithinRadius). See [`mean`](@ref mean(::AbstractManifold, ::AbstractVector, ::AbstractVector, ::GeodesicInterpolation)) for a description of `kwargs`. @@ -388,7 +388,7 @@ end ) Estimate the Riemannian center of mass of `x` using -[`ExtrinsicEstimation`](@ref), i.e. by computing the mean in the embedding and projecting +[`ExtrinsicEstimation`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/functions/#ManifoldsBase.ExtrinsicEstimation), i.e. by computing the mean in the embedding and projecting the result back. See [`mean`](@ref mean(::AbstractManifold, ::AbstractVector, ::AbstractVector, ::GeodesicInterpolation)) @@ -448,7 +448,7 @@ Compute the (optionally weighted) Riemannian median of the vector `x` of points where $\mathrm{d}_{\mathcal M}$ denotes the Riemannian [`distance`](@ref). This function is nonsmooth (i.e nondifferentiable). -In the general case, the [`CyclicProximalPointEstimation`](@ref) is used to compute the +In the general case, the [`CyclicProximalPointEstimation`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/functions/#ManifoldsBase.CyclicProximalPointEstimation) is used to compute the median. However, this default may be overloaded for specific manifolds. Compute the median using the specified `method`. @@ -479,7 +479,7 @@ end kwargs..., ) -Compute the median using [`CyclicProximalPointEstimation`](@ref). +Compute the median using [`CyclicProximalPointEstimation`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/functions/#ManifoldsBase.CyclicProximalPointEstimation). Optionally, provide `p0`, the starting point (by default set to the first data point). `stop_iter` denotes the maximal number of iterations to perform @@ -510,7 +510,7 @@ Statistics.median( kwargs..., ) -Estimate the median of `x` using [`ExtrinsicEstimation`](@ref), i.e. by computing the median +Estimate the median of `x` using [`ExtrinsicEstimation`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/functions/#ManifoldsBase.ExtrinsicEstimation), i.e. by computing the median in the embedding and projecting the result back. See [`median`](@ref median(::AbstractManifold, ::AbstractVector, ::AbstractVector, ::CyclicProximalPointEstimation)) @@ -537,7 +537,7 @@ Statistics.median( kwargs..., ) -Compute the median using [`WeiszfeldEstimation`](@ref). +Compute the median using [`WeiszfeldEstimation`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/functions/#ManifoldsBase.WeiszfeldEstimation). Optionally, provide `p0`, the starting point (by default set to the first data point). `stop_iter` denotes the maximal number of iterations to perform @@ -888,7 +888,7 @@ If `shuffle_rng` is provided, it is used to shuffle the order in which the points are considered. Optionally, pass `retraction` and `inverse_retraction` method types to specify the (inverse) retraction. -See [`GeodesicInterpolation`](@ref) for details on the geodesic +See [`GeodesicInterpolation`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/functions/#ManifoldsBase.GeodesicInterpolation) for details on the geodesic interpolation method. !!! note @@ -954,7 +954,7 @@ end Use repeated weighted geodesic interpolation to estimate the mean. Simultaneously, use a Welford-like recursion to estimate the variance. -See [`GeodesicInterpolationWithinRadius`](@ref) and +See [`GeodesicInterpolationWithinRadius`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/functions/#ManifoldsBase.GeodesicInterpolationWithinRadius) and [`mean_and_var`](@ref mean_and_var(::AbstractManifold, ::AbstractVector, ::AbstractWeights, ::GeodesicInterpolation)) for more information. """ diff --git a/tutorials/getstarted.qmd b/tutorials/getstarted.qmd index 9f47a8e470..283d79f777 100644 --- a/tutorials/getstarted.qmd +++ b/tutorials/getstarted.qmd @@ -429,7 +429,7 @@ Again, for all of these, the concrete types only have to be used if you want to ````{=commonmark} ```@bibliography -Pages = ["tutorials/getstarted.md"] +Pages = ["getstarted.md"] Canonical=false ``` ```` \ No newline at end of file diff --git a/tutorials/integration.qmd b/tutorials/integration.qmd index 019ac52199..29ae261dad 100644 --- a/tutorials/integration.qmd +++ b/tutorials/integration.qmd @@ -175,6 +175,6 @@ See [BoyaSudarshanTilma:2003](@cite) for some considerations on symmetric spaces ## Literature ```@bibliography -Pages = ["tutorials/integration.md"] +Pages = ["integration.md"] Canonical=false ``` \ No newline at end of file From e2bc5bb3217d3e3ddb61d24c32adca6adb0d11f2 Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Wed, 27 Dec 2023 09:56:55 +0100 Subject: [PATCH 3/5] Fix an edge case of PT on Circle (#695) --- NEWS.md | 131 +++++++++++++++++++++------------------ Project.toml | 2 +- src/manifolds/Circle.jl | 2 +- test/manifolds/circle.jl | 5 +- 4 files changed, 75 insertions(+), 65 deletions(-) diff --git a/NEWS.md b/NEWS.md index 123ea7737b..780bf3cc6a 100644 --- a/NEWS.md +++ b/NEWS.md @@ -5,6 +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.10] – 2023-12-27 + +### Fixed + +* Fixed mixed array index number in-place `parallel_transport_to!` on real `Circle`, on Julia 1.6. +* Violations of MD004 lint rule in this file. + ## [0.9.9] – 2023-12-25 ### Fixed @@ -18,115 +25,115 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Fixed -- Improved distribution of random vector generation for rotation matrices and complex circle. +* Improved distribution of random vector generation for rotation matrices and complex circle. ## [0.9.7] - 2023-11-14 ### Fixed -- Fixed `is_flat` for `CholeskySpace` and `SymmetricPositiveDefinite` with `LogCholeskyMetric` [https://github.com/JuliaManifolds/Manifolds.jl/issues/684](https://github.com/JuliaManifolds/Manifolds.jl/issues/684). +* Fixed `is_flat` for `CholeskySpace` and `SymmetricPositiveDefinite` with `LogCholeskyMetric` [https://github.com/JuliaManifolds/Manifolds.jl/issues/684](https://github.com/JuliaManifolds/Manifolds.jl/issues/684). ## [0.9.6] - 2023-11-09 ### Fixed -- Fixed real coefficient basis for complex circle (an issue exposed by [https://github.com/JuliaManifolds/ManifoldsBase.jl/pull/173](https://github.com/JuliaManifolds/ManifoldsBase.jl/pull/173)). -- Fixed `VeeOrthogonalBasis` test for non-real manifolds. +* Fixed real coefficient basis for complex circle (an issue exposed by [https://github.com/JuliaManifolds/ManifoldsBase.jl/pull/173](https://github.com/JuliaManifolds/ManifoldsBase.jl/pull/173)). +* Fixed `VeeOrthogonalBasis` test for non-real manifolds. ## [0.9.5] - 2023-11-08 ### Changed -- `identity_element` now returns a complex matrix for unitary group. -- `number_of_coordinates` is now exported. +* `identity_element` now returns a complex matrix for unitary group. +* `number_of_coordinates` is now exported. ## [0.9.4] - 2023-11-06 ### Added -- Functions `inv_diff`, `inv_diff!`, `adjoint_inv_diff` and `adjoint_inv_diff!` that correspond to differentials and pullbacks of group inversion. -- Julia 1.10-rc CI workflow. +* Functions `inv_diff`, `inv_diff!`, `adjoint_inv_diff` and `adjoint_inv_diff!` that correspond to differentials and pullbacks of group inversion. +* Julia 1.10-rc CI workflow. ### Changed -- Documentation project files are marked as compatible with `BoundaryValueDiffEq` v5. +* Documentation project files are marked as compatible with `BoundaryValueDiffEq` v5. ### Fixed -- Fixed issue with incorrect implementation of `apply_diff_group` in `GroupOperationAction` with left backward and right forward action [#669](https://github.com/JuliaManifolds/Manifolds.jl/issues/669). +* Fixed issue with incorrect implementation of `apply_diff_group` in `GroupOperationAction` with left backward and right forward action [#669](https://github.com/JuliaManifolds/Manifolds.jl/issues/669). ## [0.9.3] - 2023-10-28 ### Added -- Support for `BoundaryValueDiffEq` v5. +* Support for `BoundaryValueDiffEq` v5. ## [0.9.2] - 2023-10-27 ### Added -- `rand(G; vector_at=Identity(G))` now works for translation, special orthogonal and special Euclidean groups `G` (issue [#665](https://github.com/JuliaManifolds/Manifolds.jl/issues/665)). -- `get_embedding` now works for `GeneralUnitaryMultiplicationGroup`. -- Github action that checks for NEWS.md changes. +* `rand(G; vector_at=Identity(G))` now works for translation, special orthogonal and special Euclidean groups `G` (issue [#665](https://github.com/JuliaManifolds/Manifolds.jl/issues/665)). +* `get_embedding` now works for `GeneralUnitaryMultiplicationGroup`. +* Github action that checks for NEWS.md changes. ## [0.9.1] - 2023-10-25 ### Added -- a new retraction and its inverse for the fixed Rank Manifolds, the orthographic retraction. +* a new retraction and its inverse for the fixed Rank Manifolds, the orthographic retraction. ## [0.9.0] - 2023-10-24 ### Added -- Vector bundles are generalized to fiber bundles. Old `BundleFibers` functionality was reworked to better match mathematical abstractions. Fiber bundle functionality is experimental and minor changes may happen without a breaking release, with the exception of `TangentBundle` which is considered to be stable. -- `RotationTranslationAction` is introduced. +* Vector bundles are generalized to fiber bundles. Old `BundleFibers` functionality was reworked to better match mathematical abstractions. Fiber bundle functionality is experimental and minor changes may happen without a breaking release, with the exception of `TangentBundle` which is considered to be stable. +* `RotationTranslationAction` is introduced. ### Changed -- Sizes of all manifolds can now be either encoded in type or stored in a field to avoid over-specialization. +* Sizes of all manifolds can now be either encoded in type or stored in a field to avoid over-specialization. The default is set to store the size in type parameter (except for `PowerManifold` and its variants), replicating the previous behavior. For field storage, pass the `parameter=:field` keyword argument to manifold constructor. For example statically sized `CenteredMatrices{m,n}` is now `CenteredMatrices{TypeParameter{Tuple{m,n}}}`, whereas the type of special Euclidean group with field-stored size is `CenteredMatrices{Tuple{Int,Int}}`. Similar change applies to: - - `CenteredMatrices{m,n}`, - - `CholeskySpace{N}`, - - `Elliptope{N,K}`, - - `Euclidean`, - - `FixedRankMatrices{m,n,k}`, - - `KendallsPreShapeSpace{n,k}`, - - `KendallsShapeSpace{n,k}`, - - `GeneralLinear{n}`, - - `GeneralUnitaryMultiplicationGroup{n}`, - - `GeneralizedGrassmann{n,k}`, - - `GeneralizedStiefel{n,k}`, - - `Grassmann{n,k}`, - - `Heisenberg{n}`, - - `Hyperbolic{n}`, - - `MultinomialMatrices{N,M}`, - - `MultinomialDoublyStochastic{n}`, - - `MultinomialSymmetric{n}`, - - `Orthogonal{n}`, - - `PowerManifold`, - - `PositiveArrays`, - - `PositiveMatrices`, - - `PositiveNumbers`, - - `ProbabilitySimplex{n}`, - - `SPDFixedDeterminant{n}`, - - `SpecialLinear{n}`, - - `SpecialOrthogonal{n}`, - - `SpecialUnitary{n}`, - - `SpecialEuclidean{n}`, - - `SpecialEuclideanManifold{n}`, - - `Spectrahedron{n,k}`, - - `SphereSymmetricMatrices{N}`, - - `Stiefel{n,k}`, - - `SymmetricMatrices{N}`, - - `SymmetricPositiveDefinite{n}`, - - `SymmetricPositiveSemidefiniteFixedRank{n,k}`, - - `Symplectic{n}`, - - `SymplecticStiefel{n,k}`, - - `TranslationGroup`, - - `Tucker`. + * `CenteredMatrices{m,n}`, + * `CholeskySpace{N}`, + * `Elliptope{N,K}`, + * `Euclidean`, + * `FixedRankMatrices{m,n,k}`, + * `KendallsPreShapeSpace{n,k}`, + * `KendallsShapeSpace{n,k}`, + * `GeneralLinear{n}`, + * `GeneralUnitaryMultiplicationGroup{n}`, + * `GeneralizedGrassmann{n,k}`, + * `GeneralizedStiefel{n,k}`, + * `Grassmann{n,k}`, + * `Heisenberg{n}`, + * `Hyperbolic{n}`, + * `MultinomialMatrices{N,M}`, + * `MultinomialDoublyStochastic{n}`, + * `MultinomialSymmetric{n}`, + * `Orthogonal{n}`, + * `PowerManifold`, + * `PositiveArrays`, + * `PositiveMatrices`, + * `PositiveNumbers`, + * `ProbabilitySimplex{n}`, + * `SPDFixedDeterminant{n}`, + * `SpecialLinear{n}`, + * `SpecialOrthogonal{n}`, + * `SpecialUnitary{n}`, + * `SpecialEuclidean{n}`, + * `SpecialEuclideanManifold{n}`, + * `Spectrahedron{n,k}`, + * `SphereSymmetricMatrices{N}`, + * `Stiefel{n,k}`, + * `SymmetricMatrices{N}`, + * `SymmetricPositiveDefinite{n}`, + * `SymmetricPositiveSemidefiniteFixedRank{n,k}`, + * `Symplectic{n}`, + * `SymplecticStiefel{n,k}`, + * `TranslationGroup`, + * `Tucker`. For example @@ -166,12 +173,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 end ``` -- Argument order for type aliases `RotationActionOnVector` and `RotationTranslationActionOnVector`: most often dispatched on argument is now first. -- A more consistent handling of action direction was introduced. 4-valued `ActionDirection` was split into 2-valued `ActionDirection` (either left or right action) and `GroupActionSide` (action acting from the left or right side). See [https://github.com/JuliaManifolds/Manifolds.jl/issues/637](https://github.com/JuliaManifolds/Manifolds.jl/issues/637) for a design discussion. +* Argument order for type aliases `RotationActionOnVector` and `RotationTranslationActionOnVector`: most often dispatched on argument is now first. +* A more consistent handling of action direction was introduced. 4-valued `ActionDirection` was split into 2-valued `ActionDirection` (either left or right action) and `GroupActionSide` (action acting from the left or right side). See [https://github.com/JuliaManifolds/Manifolds.jl/issues/637](https://github.com/JuliaManifolds/Manifolds.jl/issues/637) for a design discussion. ### Removed -- `ProductRepr` is removed; please use `ArrayPartition` instead. -- Default methods throwing "not implemented" `ErrorException` for some group-related operations. Standard `MethodError` is now thrown instead. -- `LinearAffineMetric` was deprecated in a previous release and the symbol is now removed. +* `ProductRepr` is removed; please use `ArrayPartition` instead. +* Default methods throwing "not implemented" `ErrorException` for some group-related operations. Standard `MethodError` is now thrown instead. +* `LinearAffineMetric` was deprecated in a previous release and the symbol is now removed. Please use `AffineInvariantMetric` instead. diff --git a/Project.toml b/Project.toml index 54bb294f41..b297892cc5 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.9" +version = "0.9.10" [deps] Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" diff --git a/src/manifolds/Circle.jl b/src/manifolds/Circle.jl index c6b2d7a2f3..a82cff544c 100644 --- a/src/manifolds/Circle.jl +++ b/src/manifolds/Circle.jl @@ -532,7 +532,7 @@ function parallel_transport_to(M::Circle{ℂ}, p::Number, X::Number, q::Number) return Y end -parallel_transport_to!(::Circle{ℝ}, Y, p, X, q) = (Y .= X) +parallel_transport_to!(::Circle{ℝ}, Y, p, X, q) = (Y .= X[]) function parallel_transport_to!(M::Circle{ℂ}, Y, p, X, q) X_pq = log(M, p, q) Xnorm = norm(M, p, X_pq) diff --git a/test/manifolds/circle.jl b/test/manifolds/circle.jl index 46bfe9d226..bc829fb91f 100644 --- a/test/manifolds/circle.jl +++ b/test/manifolds/circle.jl @@ -288,10 +288,13 @@ using Manifolds: TFVector, CoTFVector ) end end - @testset "Mixed array dimensions for exp" begin + @testset "Mixed array dimensions for exp and PT" begin + # this is an issue on Julia 1.6 but not later releases M = Circle() p = fill(0.0) exp!(M, p, p, [1.0], 2.0) @test p ≈ fill(2.0) + parallel_transport_to!(M, p, p, [4.0], p) + @test p ≈ fill(4.0) end end From fab7b4bd1336a4ea2d283ca9041f8e7b3faca0c4 Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Wed, 27 Dec 2023 13:00:30 +0100 Subject: [PATCH 4/5] Compatibility with RecursiveArrayTools v3 (#696) --- NEWS.md | 4 ++++ Project.toml | 2 +- test/manifolds/product_manifold.jl | 4 ++-- 3 files changed, 7 insertions(+), 3 deletions(-) diff --git a/NEWS.md b/NEWS.md index 780bf3cc6a..16e457f554 100644 --- a/NEWS.md +++ b/NEWS.md @@ -7,6 +7,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [0.9.10] – 2023-12-27 +### Added + +* Compatibility with `RecursiveArrayTools` v3. + ### Fixed * Fixed mixed array index number in-place `parallel_transport_to!` on real `Circle`, on Julia 1.6. diff --git a/Project.toml b/Project.toml index b297892cc5..0574f8379e 100644 --- a/Project.toml +++ b/Project.toml @@ -59,7 +59,7 @@ Plots = "1" Quaternions = "0.5, 0.6, 0.7" Random = "1.6" RecipesBase = "1.1" -RecursiveArrayTools = "2" +RecursiveArrayTools = "2, 3" Requires = "0.5, 1" SimpleWeightedGraphs = "1.2" SpecialFunctions = "0.8, 0.9, 0.10, 1.0, 2" diff --git a/test/manifolds/product_manifold.jl b/test/manifolds/product_manifold.jl index f9233a984e..5bbf9cd3cf 100644 --- a/test/manifolds/product_manifold.jl +++ b/test/manifolds/product_manifold.jl @@ -350,12 +350,12 @@ using RecursiveArrayTools: ArrayPartition Z = -X H1 = riemannian_Hessian(N, p, Y, Z, X) H2 = ArrayPartition( - [riemannian_Hessian(M, p[i, :], Y[i, :], Z[i, :], X[i, :]) for i in 1:2]..., + [riemannian_Hessian(M, p.x[i], Y.x[i], Z.x[i], X.x[i]) for i in 1:2]..., ) @test H1 == H2 V = ArrayPartition([0.2, 0.0, 0.0], [0.0, 0.0, 0.3]) W1 = Weingarten(N, p, X, V) - W2 = ArrayPartition([Weingarten(M, p[i, :], X[i, :], V[i, :]) for i in 1:2]...) + W2 = ArrayPartition([Weingarten(M, p.x[i], X.x[i], V.x[i]) for i in 1:2]...) @test W1 == W2 end @testset "Manifold volume" begin From 81a43d434a1ed6f515c00dd9814fd45d0d99ae2e Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Wed, 27 Dec 2023 18:37:16 +0100 Subject: [PATCH 5/5] Fix 0-index Euclidean just like Circle (#697) --- NEWS.md | 6 ++++++ Project.toml | 2 +- src/manifolds/Euclidean.jl | 6 ++++++ test/manifolds/euclidean.jl | 11 +++++++++++ 4 files changed, 24 insertions(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index 16e457f554..0167ac539b 100644 --- a/NEWS.md +++ b/NEWS.md @@ -5,6 +5,12 @@ 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.11] – 2023-12-27 + +### Fixed + +* Fixed mixed array index number in-place `parallel_transport_to!` on zero-index `Euclidean`, on Julia 1.6. + ## [0.9.10] – 2023-12-27 ### Added diff --git a/Project.toml b/Project.toml index 0574f8379e..0f60e14baf 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.10" +version = "0.9.11" [deps] Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" diff --git a/src/manifolds/Euclidean.jl b/src/manifolds/Euclidean.jl index f2160ef492..4dc39a7da6 100644 --- a/src/manifolds/Euclidean.jl +++ b/src/manifolds/Euclidean.jl @@ -213,6 +213,8 @@ Base.exp(::Euclidean, p, X, t::Number) = p .+ t .* X exp!(::Euclidean, q, p, X) = (q .= p .+ X) exp!(::Euclidean, q, p, X, t::Number) = (q .= p .+ t .* X) +exp!(::Euclidean{TypeParameter{Tuple{}}}, q, p, X, t::Number) = (q .= p[] + t * X[]) +exp!(::Euclidean{Tuple{}}, q, p, X, t::Number) = (q .= p[] + t * X[]) function get_basis_diagonalizing( M::Euclidean, @@ -677,6 +679,10 @@ the parallel transport on [`Euclidean`](@ref) is the identiy, i.e. returns `X`. """ parallel_transport_to(::Euclidean, ::Any, X, ::Any) = X parallel_transport_to!(::Euclidean, Y, ::Any, X, ::Any) = copyto!(Y, X) +function parallel_transport_to!(::Euclidean{TypeParameter{Tuple{}}}, Y, ::Any, X, ::Any) + return copyto!(Y, X[]) +end +parallel_transport_to!(::Euclidean{Tuple{}}, Y, ::Any, X, ::Any) = copyto!(Y, X[]) @doc raw""" project(M::Euclidean, p) diff --git a/test/manifolds/euclidean.jl b/test/manifolds/euclidean.jl index c92fe16ce4..e83572c67e 100644 --- a/test/manifolds/euclidean.jl +++ b/test/manifolds/euclidean.jl @@ -446,4 +446,15 @@ using FiniteDifferences ManifoldDiff.βdifferential_shortest_geodesic_startpoint, ) === 2.0 end + + @testset "Mixed array dimensions for exp and PT" begin + # this is an issue on Julia 1.6 but not later releases + for M in [Euclidean(), Euclidean(; parameter=:field)] + p = fill(0.0) + exp!(M, p, p, [1.0], 2.0) + @test p ≈ fill(2.0) + parallel_transport_to!(M, p, p, [4.0], p) + @test p ≈ fill(4.0) + end + end end