Improve naming, test coverage and REPL display.

Author: kellertuer

Changelog.md

@@ -12,7 +12,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 * a `ScaledManifoldObjective` to easier build scaled versions of objectives,
   especially turn maximisation problems into minimisation ones using a scaling of `-1`.
-* Introduce a `ConstrainedSetObjective`
+* Introduce a `ManifoldConstrainedSetObjective`
 * Introduce a `projected_gradient_method`
 
 

docs/src/plans/objective.md

@@ -205,7 +205,7 @@ linearized_forward_operator
 
 ```@docs
 ConstrainedManifoldObjective
-ConstrainedSetObjective
+ManifoldConstrainedSetObjective
 ```
 
 It might be beneficial to use the adapted problem to specify different ranges for the gradients of the constraints

src/Manopt.jl

@@ -297,6 +297,7 @@ export AbstractDecoratedManifoldObjective,
     AbstractManifoldSubObjective,
     AbstractPrimalDualManifoldObjective,
     ConstrainedManifoldObjective,
+    ManifoldConstrainedSetObjective,
     EmbeddedManifoldObjective,
     ScaledManifoldObjective,
     ManifoldCountObjective,

src/plans/constrained_set_plan.jl

@@ -1,5 +1,5 @@
 """
-    ConstrainedSetObjective{E, MO, PF, IF} <: AbstractManifoldObjective{E}
+    ManifoldConstrainedSetObjective{E, MO, PF, IF} <: AbstractManifoldObjective{E}
 
 Model a constrained objective restricted to a set
 
@@ -18,7 +18,7 @@ where ``$(_tex(:Cal,"C")) ⊂ $(_math(:M))`` is a convex closed subset.
 
 # Constructor
 
-    ConstrainedSetObjective(f, grad_f, project!!; kwargs...)
+    ManifoldConstrainedSetObjective(f, grad_f, project!!; kwargs...)
 
 Generate the constrained objective for a given function `f` its gradient `grad_f` and a projection `project!!` ``$(_tex(:proj))_{$(_tex(:Cal,"C"))}``.
 
@@ -28,22 +28,24 @@ $(_var(:Keyword, :evaluation))
 * `indicator=nothing`: the indicator function ``ι_{$(_tex(:Cal,"C"))}(p)``. If not provided a test, whether the projection yields the same point is performed.
   For the [`InplaceEvaluation`](@ref) this required one allocation.
 """
-struct ConstrainedSetObjective{
+struct ManifoldConstrainedSetObjective{
     E<:AbstractEvaluationType,MO<:AbstractManifoldObjective,PF,IF
 } <: AbstractManifoldObjective{E}
     objective::MO
     project!!::PF
     indicator::IF
 end
 
-function ConstrainedSetObjective(
+function ManifoldConstrainedSetObjective(
     f, grad_f, project!!::PF; evaluation::E=AllocatingEvaluation(), indicator=nothing
 ) where {PF,E<:AbstractEvaluationType}
     obj = ManifoldGradientObjective(f, grad_f; evaluation=evaluation)
     if isnothing(indicator)
         if evaluation isa AllocatingEvaluation
             ind(M, p) = (distance(M, p, project!!(M, p)) ≈ 0 ? 0 : Inf)
-            return ConstrainedSetObjective{E,typeof(obj),typeof(project!!),typeof(ind)}(
+            return ManifoldConstrainedSetObjective{
+                E,typeof(obj),typeof(project!!),typeof(ind)
+            }(
                 obj, project!!, ind
             )
         elseif evaluation isa InplaceEvaluation
@@ -52,39 +54,43 @@ function ConstrainedSetObjective(
                 project!!(M, q, p)
                 return distance(M, p, q) ≈ 0 ? 0 : Inf
             end
-            return ConstrainedSetObjective{E,typeof(obj),typeof(project!!),typeof(ind)}(
+            return ManifoldConstrainedSetObjective{
+                E,typeof(obj),typeof(project!!),typeof(ind)
+            }(
                 obj, project!!, ind
             )
         end
     end
-    return ConstrainedSetObjective{E,typeof(obj),typeof(project!!),typeof(indicator)}(
+    return ManifoldConstrainedSetObjective{
+        E,typeof(obj),typeof(project!!),typeof(indicator)
+    }(
         obj, project!!, indicator
     )
 end
 
-function get_cost(M::AbstractManifold, cso::ConstrainedSetObjective, p)
+function get_cost(M::AbstractManifold, cso::ManifoldConstrainedSetObjective, p)
     return get_cost(M, cso.objective, p)
 end
-function get_cost_function(cso::ConstrainedSetObjective, recursive=false)
-    return get_cost_function(cso.objective)
+function get_cost_function(cso::ManifoldConstrainedSetObjective, recursive=false)
+    return get_cost_function(cso.objective, recursive)
 end
-function get_gradient_function(cso::ConstrainedSetObjective, recursive=false)
-    return get_gradient_function(cso.objective)
+function get_gradient_function(cso::ManifoldConstrainedSetObjective, recursive=false)
+    return get_gradient_function(cso.objective, recursive)
 end
-function get_gradient(M::AbstractManifold, cso::ConstrainedSetObjective, p)
+function get_gradient(M::AbstractManifold, cso::ManifoldConstrainedSetObjective, p)
     return get_gradient(M, cso.objective, p)
 end
-function get_gradient!(M::AbstractManifold, X, cso::ConstrainedSetObjective, p)
+function get_gradient!(M::AbstractManifold, X, cso::ManifoldConstrainedSetObjective, p)
     return get_gradient!(M, X, cso.objective, p)
 end
 
 _doc_get_projected_point = """
     get_projected_point(amp::AbstractManoptProblem, p)
     get_projected_point!(amp::AbstractManoptProblem, q, p)
-    get_projected_point(M::AbstractManifold, cso::ConstrainedSetObjective, p)
-    get_projected_point!(M::AbstractManifold, q, cso::ConstrainedSetObjective, p)
+    get_projected_point(M::AbstractManifold, cso::ManifoldConstrainedSetObjective, p)
+    get_projected_point!(M::AbstractManifold, q, cso::ManifoldConstrainedSetObjective, p)
 
-Project `p` with the projection that is stored within the [`ConstrainedSetObjective`](@ref).
+Project `p` with the projection that is stored within the [`ManifoldConstrainedSetObjective`](@ref).
 This can be done in-place of `q`.
 """
 
@@ -98,29 +104,29 @@ function get_projected_point!(amp::AbstractManoptProblem, q, p)
 end
 
 @doc "$(_doc_get_projected_point)"
-get_projected_point(M::AbstractManifold, cso::ConstrainedSetObjective, p)
+get_projected_point(M::AbstractManifold, cso::ManifoldConstrainedSetObjective, p)
 function get_projected_point(
-    M::AbstractManifold, cso::ConstrainedSetObjective{AllocatingEvaluation}, p
+    M::AbstractManifold, cso::ManifoldConstrainedSetObjective{AllocatingEvaluation}, p
 )
     return cso.project!!(M, p)
 end
 function get_projected_point(
-    M::AbstractManifold, cso::ConstrainedSetObjective{InplaceEvaluation}, p
+    M::AbstractManifold, cso::ManifoldConstrainedSetObjective{InplaceEvaluation}, p
 )
     q = copy(M, p)
     cso.project!!(M, q, p)
     return q
 end
 @doc "$(_doc_get_projected_point)"
-get_projected_point!(M::AbstractManifold, q, cso::ConstrainedSetObjective, p)
+get_projected_point!(M::AbstractManifold, q, cso::ManifoldConstrainedSetObjective, p)
 function get_projected_point!(
-    M::AbstractManifold, q, cso::ConstrainedSetObjective{AllocatingEvaluation}, p
+    M::AbstractManifold, q, cso::ManifoldConstrainedSetObjective{AllocatingEvaluation}, p
 )
     copyto!(M, q, cso.project!!(M, p))
     return q
 end
 function get_projected_point!(
-    M::AbstractManifold, q, cso::ConstrainedSetObjective{InplaceEvaluation}, p
+    M::AbstractManifold, q, cso::ManifoldConstrainedSetObjective{InplaceEvaluation}, p
 )
     cso.project!!(M, q, p)
     return q

src/plans/gradient_plan.jl

@@ -92,8 +92,10 @@ function ManifoldCostGradientObjective(
     return ManifoldCostGradientObjective{typeof(evaluation),CG}(costgrad)
 end
 
-get_cost_function(cgo::ManifoldCostGradientObjective) = (M, p) -> get_cost(M, cgo, p)
-function get_gradient_function(cgo::ManifoldCostGradientObjective)
+function get_cost_function(cgo::ManifoldCostGradientObjective, recursive=false)
+    return (M, p) -> get_cost(M, cgo, p)
+end
+function get_gradient_function(cgo::ManifoldCostGradientObjective, recursive=false)
     return (M, p) -> get_gradient(M, cgo, p)
 end
 

src/solvers/projected_gradient_method.jl

@@ -142,20 +142,25 @@ function get_reason(c::StopWhenProjectedGradientStationary)
     end
     return ""
 end
+indicates_convergence(c::StopWhenProjectedGradientStationary) = true
+function show(io::IO, c::StopWhenProjectedGradientStationary)
+    return print(
+        io, "StopWhenProjectedGradientStationary($(c.threshold))\n    $(status_summary(c))"
+    )
+end
 function status_summary(c::StopWhenProjectedGradientStationary)
     has_stopped = (c.at_iteration >= 0)
     s = has_stopped ? "reached" : "not reached"
-    return "stopped after $(c.threshold):\t$s"
+    return "projected gradient stationary (<$(c.threshold)): \t$s"
 end
-
 #
 #
 # The solver
 _doc_pgm = """
     projected_gradient_method(M, f, grad_f, proj, p=rand(M); kwargs...)
-    projected_gradient_method(M, obj::ConstrainedSetObjective, p=rand(M); kwargs...)
+    projected_gradient_method(M, obj::ManifoldConstrainedSetObjective, p=rand(M); kwargs...)
     projected_gradient_method!(M, f, grad_f, proj, p; kwargs...)
-    projected_gradient_method!(M, obj::ConstrainedSetObjective, p; kwargs...)
+    projected_gradient_method!(M, obj::ManifoldConstrainedSetObjective, p; kwargs...)
 
 Compute the projected gradient method for the constrained problem
 
@@ -202,12 +207,12 @@ end
 function projected_gradient_method(
     M, f, grad_f, proj, p; indicator=nothing, evaluation=AllocatingEvaluation(), kwargs...
 )
-    cs_obj = ConstrainedSetObjective(
+    cs_obj = ManifoldConstrainedSetObjective(
         f, grad_f, proj; evaluation=evaluation, indicator=indicator
     )
     return projected_gradient_method(M, cs_obj, p; kwargs...)
 end
-function projected_gradient_method(M, obj::ConstrainedSetObjective, p; kwargs...)
+function projected_gradient_method(M, obj::ManifoldConstrainedSetObjective, p; kwargs...)
     q = copy(M, p)
     return projected_gradient_method!(M, obj, q; kwargs...)
 end
@@ -216,14 +221,14 @@ end
 function projected_gradient_method!(
     M, f, grad_f, proj, p; indicator=nothing, evaluation=AllocatingEvaluation(), kwargs...
 )
-    cs_obj = ConstrainedSetObjective(
+    cs_obj = ManifoldConstrainedSetObjective(
         f, grad_f, proj; evaluation=evaluation, indicator=indicator
     )
     return projected_gradient_method!(M, cs_obj, p; kwargs...)
 end
 function projected_gradient_method!(
     M,
-    obj::ConstrainedSetObjective,
+    obj::ManifoldConstrainedSetObjective,
     p;
     backtrack::Stepsize=ArmijoLinesearchStepsize(M; stop_increasing_at_step=0),
     retraction_method::AbstractRetractionMethod=default_retraction_method(M, typeof(p)),

test/plans/test_constrained_set_plan.jl

@@ -0,0 +1,88 @@
+using Manifolds, Manopt, Random, Test
+
+@testset "Constained set objective" begin
+    M = Hyperbolic(2)
+    c = Manifolds._hyperbolize(M, [0, 0])
+    r = 1.0
+    N = 200
+    σ = 1.5
+    Random.seed!(42)
+    # N random points moved to top left to have a mean outside
+    pts = [
+        exp(
+            M,
+            c,
+            get_vector(
+                M,
+                c,
+                σ .* randn(manifold_dimension(M)) .+ [2.5, 2.5],
+                DefaultOrthonormalBasis(),
+            ),
+        ) for _ in 1:N
+    ]
+    f(M, p) = 1 / (2 * length(pts)) .* sum(distance(M, p, q)^2 for q in pts)
+    grad_f(M, p) = -1 / length(pts) .* sum(log(M, p, q) for q in pts)
+    function grad_f!(M, X, p)
+        zero_vector!(M, X, p)
+        Y = zero_vector(M, p)
+        for q in pts
+            log!(M, Y, p, q)
+            X .+= Y
+        end
+        X .*= -1 / length(pts)
+        return X
+    end
+    function project_C(M, p)
+        X = log(M, c, p)
+        n = norm(M, c, X)
+        q = (n > r) ? exp(M, c, (r / n) * X) : copy(M, p)
+        return q
+    end
+    function project_C!(M, q, p; X=zero_vector(M, c))
+        log!(M, X, c, p)
+        n = norm(M, c, X)
+        if (n > r)
+            exp!(M, q, c, (r / n) * X)
+        else
+            copyto!(M, q, p)
+        end
+        return q
+    end
+    g(M, p) = distance(M, c, p)^2 - r^2
+    indicator_C(M, p) = (g(M, p) ≤ 0) ? 0 : Inf
+
+    csoa = ManifoldConstrainedSetObjective(f, grad_f, project_C)
+    csoa2 = ManifoldConstrainedSetObjective(f, grad_f, project_C; indicator=indicator_C)
+    csoi = ManifoldConstrainedSetObjective(
+        f, grad_f!, project_C!; evaluation=InplaceEvaluation()
+    )
+    csoi2 = ManifoldConstrainedSetObjective(
+        f, grad_f!, project_C!; evaluation=InplaceEvaluation(), indicator=indicator_C
+    )
+
+    for objective in [csoa, csoa2, csoi, csoi2]
+        @test get_cost(M, objective, c) == f(M, c)
+        @test Manopt.get_cost_function(objective)(M, c) == f(M, c)
+        @test get_gradient(M, objective, c) == grad_f(M, c)
+        X = zero_vector(M, c)
+        get_gradient!(M, X, objective, c)
+        Y = zero_vector(M, c)
+        grad_f!(M, Y, c)
+        @test X == Y
+        if objective ∈ [csoa, csoa2]
+            @test Manopt.get_gradient_function(objective)(M, c) == grad_f(M, c)
+        else
+            Manopt.get_gradient_function(objective)(M, X, c) == grad_f!(M, Y, c)
+            @test X == Y
+        end
+        dmp = DefaultManoptProblem(M, objective)
+        p = get_projected_point(dmp, c)
+        @test p == c # c is already in C
+        get_projected_point!(dmp, p, c)
+        @test p == c
+        p = get_projected_point(M, objective, c)
+        @test p == c
+        get_projected_point!(M, p, objective, c)
+        @test p == c
+    end
+end

test/runtests.jl

@@ -20,6 +20,7 @@ include("utils/example_tasks.jl")
         include("plans/test_nonlinear_least_squares_plan.jl")
         include("plans/test_gradient_plan.jl")
         include("plans/test_constrained_plan.jl")
+        include("plans/test_constrained_set_plan.jl")
         include("plans/test_hessian_plan.jl")
         include("plans/test_parameters.jl")
         include("plans/test_primal_dual_plan.jl")

test/solvers/test_projected_gradient.jl

@@ -57,16 +57,35 @@ using Manifolds, Manopt, Random, Test
         stopping_criterion=StopAfterIteration(150) |
                            StopWhenProjectedGradientStationary(M, 1e-7),
     )
-    mean_pg_2 = copy(M, c)
-    projected_gradient_method!(
+    Random.seed!(42)
+    mean_pg_2 = projected_gradient_method(
+        M,
+        f,
+        grad_f,
+        project_C;
+        stopping_criterion=StopAfterIteration(150) |
+                           StopWhenProjectedGradientStationary(M, 1e-7),
+    )
+    @test isapprox(M, mean_pg_1, mean_pg_2)
+    mean_pg_3 = copy(M, c)
+    st = projected_gradient_method!(
         M,
         f,
         grad_f!,
         project_C!,
-        mean_pg_2;
+        mean_pg_3;
         evaluation=InplaceEvaluation(),
         stopping_criterion=StopAfterIteration(150) |
                            StopWhenProjectedGradientStationary(M, 1e-7),
+        return_state=true,
     )
-    @test isapprox(M, mean_pg_1, mean_pg_2)
+    @test isapprox(M, mean_pg_1, mean_pg_3)
+    @test startswith(
+        repr(st), "# Solver state for `Manopt.jl`s Projected Gradient Method\n"
+    )
+    stop_when_stationary = st.stop.criteria[2]
+    @test repr(stop_when_stationary) ==
+        "StopWhenProjectedGradientStationary($(stop_when_stationary.threshold))\n    $(Manopt.status_summary(
+        stop_when_stationary))"
+    @test length(get_reason(stop_when_stationary)) > 0
 end