Introduce a :Subsolver debug symbol and a DebugWhenActive

Project.toml

@@ -1,7 +1,7 @@
 name = "Manopt"
 uuid = "0fc0a36d-df90-57f3-8f93-d78a9fc72bb5"
 authors = ["Ronny Bergmann <manopt@ronnybergmann.net>"]
-version = "0.4.30"
+version = "0.4.31"
 
 [deps]
 ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4"

docs/make.jl

@@ -86,6 +86,7 @@ makedocs(
             "Speedup using Inplace computations" => "tutorials/InplaceGradient.md",
             "Use Automatic Differentiation" => "tutorials/AutomaticDifferentiation.md",
             "Count and use a Cache" => "tutorials/CountAndCache.md",
+            "Perform Debug Output" => "tutorials/HowToDebug.md",
             "Record values" => "tutorials/HowToRecord.md",
             "Implement a Solver" => "tutorials/ImplementASolver.md",
             "Do Contrained Optimization" => "tutorials/ConstrainedOptimization.md",

docs/src/plans/state.md

@@ -98,3 +98,10 @@ AbstractGradientSolverState
 AbstractHessianSolverState
 AbstractPrimalDualSolverState
 ```
+
+For the sub problem state, there are two access functions
+
+```@docs
+get_sub_problem
+get_sub_state
+```

docs/src/tutorials/HowToDebug.md

@@ -0,0 +1,173 @@
+# How to Print Debug Output
+Ronny Bergmann
+
+This tutorial aims to illustrate how to perform debug output. For that we consider an
+example that includes a subsolver, to also consider their debug capabilities.
+
+The problem itself is hence not the main focus.
+
+We consider a nonnegative PCA which we can write as a constraint problem on the Sphere
+
+Let’s first load the necessary packages.
+
+``` julia
+using Manopt, Manifolds, Random, LinearAlgebra
+Random.seed!(42);
+```
+
+``` julia
+d = 4
+M = Sphere(d - 1)
+v0 = project(M, [ones(2)..., zeros(d - 2)...])
+Z = v0 * v0'
+#Cost and gradient
+f(M, p) = -tr(transpose(p) * Z * p) / 2
+grad_f(M, p) = project(M, p, -transpose.(Z) * p / 2 - Z * p / 2)
+# Constraints
+g(M, p) = -p # i.e. p ≥ 0
+mI = -Matrix{Float64}(I, d, d)
+# Vector of gradients of the constraint components
+grad_g(M, p) = [project(M, p, mI[:, i]) for i in 1:d]
+```
+
+Then we can take a starting point
+
+``` julia
+p0 = project(M, [ones(2)..., zeros(d - 3)..., 0.1])
+```
+
+## Simple debug output
+
+Any solver accepts the keyword `debug=`, which in the simplest case can be set to an array of strings, symbols and a number.
+
+- Strings are printed in every iteration as is (cf. [`DebugDivider`](@ref)) and should be used to finish the array with a line break.
+- the last number in the array is used with [`DebugEvery`](@ref) to print the debug only every $i$th iteration.
+- Any Symbol is converted into certain debug prints
+
+Certain symbols starting with a capital letter are mapped to certain prints, e.g. `:Cost` is mapped to [`DebugCost`](@ref)`()` to print the current cost function value. A full list is provided in the [`DebugActionFactory`](@ref).
+A special keyword is `:Stop`, which is only added to the final debug hook to print the stopping criterion.
+
+Any symbol with a small letter is mapped to fields of the [`AbstractManoptSolverState`](@ref) which is used. This way you can easily print internal data, if you know their names.
+
+Let’s look at an example first: If we want to print the current iteration number, the current cost function value as well as the value `ϵ` from the [`ExactPenaltyMethodState`](@ref). To keep the amount of print at a reasonable level, we want to only print the debug every 25th iteration.
+
+Then we can write
+
+``` julia
+p1 = exact_penalty_method(
+    M, f, grad_f, p0; g=g, grad_g=grad_g,
+    debug = [:Iteration, :Cost, " | ", :ϵ, 25, "\n", :Stop]
+);
+```
+
+    Initial f(x): -0.497512 | ϵ: 0.001
+    # 25    f(x): -0.499449 | ϵ: 0.0001778279410038921
+    # 50    f(x): -0.499995 | ϵ: 3.1622776601683734e-5
+    # 75    f(x): -0.500000 | ϵ: 5.623413251903474e-6
+    # 100   f(x): -0.500000 | ϵ: 1.0e-6
+    The value of the variable (ϵ) is smaller than or equal to its threshold (1.0e-6).
+    The algorithm performed a step with a change (6.5347623783315016e-9) less than 1.0e-6.
+
+## Advanced Debug output
+
+There is two more advanced variants that can be used. The first is a tuple of a symbol and a string, where the string is used as the format print, that most [`DebugAction`](@ref)s have. The second is, to directly provide a `DebugAction`.
+
+We can for example change the way the `:ϵ` is printed by adding a format string
+and use [`DebugCost`](@ref)`()` which is equivalent to using `:Cost`.
+Especially with the format change, the lines are more coniststent in length.
+
+``` julia
+p2 = exact_penalty_method(
+    M, f, grad_f, p0; g=g, grad_g=grad_g,
+    debug = [:Iteration, DebugCost(), (:ϵ," | ϵ: %.8f"), 25, "\n", :Stop]
+);
+```
+
+    Initial f(x): -0.497512 | ϵ: 0.00100000
+    # 25    f(x): -0.499449 | ϵ: 0.00017783
+    # 50    f(x): -0.499995 | ϵ: 0.00003162
+    # 75    f(x): -0.500000 | ϵ: 0.00000562
+    # 100   f(x): -0.500000 | ϵ: 0.00000100
+    The value of the variable (ϵ) is smaller than or equal to its threshold (1.0e-6).
+    The algorithm performed a step with a change (6.5347623783315016e-9) less than 1.0e-6.
+
+You can also write your own [`DebugAction`](@ref) functor, where the function to implement has the same signature as the `step` function, that is an [`AbstractManoptProblem`](@ref), an [`AbstractManoptSolverState`](@ref), as well as the current iterate. For example the already mentioned \[`DebugDivider](@ref)`(s)\` is given as
+
+``` julia
+mutable struct DebugDivider{TIO<:IO} <: DebugAction
+    io::TIO
+    divider::String
+    DebugDivider(divider=" | "; io::IO=stdout) = new{typeof(io)}(io, divider)
+end
+function (d::DebugDivider)(::AbstractManoptProblem, ::AbstractManoptSolverState, i::Int)
+    (i >= 0) && (!isempty(d.divider)) && (print(d.io, d.divider))
+    return nothing
+end
+```
+
+or you could implement that of course just for your specific problem or state.
+
+## Subsolver Debug
+
+most subsolvers have a `sub_kwargs` keyword, such that you can pass keywords to the sub solver as well. This works well if you do not plan to change the subsolver. If you do you can wrap your own `solver_state=` argument in a [`decorate_state!`](@ref) and pass a `debug=` password to this function call.
+Keywords in a keyword have to be passed as pairs (`:debug => [...]`).
+
+A main problem now is, that this debug is issued every sub solver call or initialisation, as the following print of just a `.` per sub solver test/call illustrates
+
+``` julia
+p3 = exact_penalty_method(
+    M, f, grad_f, p0; g=g, grad_g=grad_g,
+    debug = ["\n",:Iteration, DebugCost(), (:ϵ," | ϵ: %.8f"), 25, "\n", :Stop],
+    sub_kwargs = [:debug => ["."]]
+);
+```
+
+
+    Initial f(x): -0.497512 | ϵ: 0.00100000
+    ........................................................
+    # 25    f(x): -0.499449 | ϵ: 0.00017783
+    ..................................................
+    # 50    f(x): -0.499995 | ϵ: 0.00003162
+    ..................................................
+    # 75    f(x): -0.500000 | ϵ: 0.00000562
+    ..................................................
+    # 100   f(x): -0.500000 | ϵ: 0.00000100
+    ....The value of the variable (ϵ) is smaller than or equal to its threshold (1.0e-6).
+    The algorithm performed a step with a change (6.5347623783315016e-9) less than 1.0e-6.
+
+The different lengths of the dotted lines come from the fact that —at least in the beginning— the subsolver performs a few steps.
+
+For this issue, there is the next symbol (similar to the `:Stop`) to indicate that a debug set is a subsolver set `:Subsolver`, which introduces a [`DebugWhenActive`](@ref) that is only activated when the outer debug is actually active, i.e. [`DebugEvery`](@ref) is active itself.
+Let’s
+
+``` julia
+p4 = exact_penalty_method(
+  M, f, grad_f, p0; g=g, grad_g=grad_g,
+  debug = [:Iteration, DebugCost(), (:ϵ," | ϵ: %.8f"), 25, "\n", :Stop],
+  sub_kwargs = [
+    :debug => [" | ", :Iteration, :Cost, "\n",:Stop, :Subsolver]
+  ]
+);
+```
+
+    Initial f(x): -0.497512 | ϵ: 0.00100000
+     | Initial f(x): -0.499127
+     | # 1     f(x): -0.499147
+    The algorithm reached approximately critical point after 1 iterations; the gradient norm (0.0002121889852717264) is less than 0.001.
+    # 25    f(x): -0.499449 | ϵ: 0.00017783
+     | Initial f(x): -0.499993
+     | # 1     f(x): -0.499994
+    The algorithm reached approximately critical point after 1 iterations; the gradient norm (1.6025009584516425e-5) is less than 0.001.
+    # 50    f(x): -0.499995 | ϵ: 0.00003162
+     | Initial f(x): -0.500000
+     | # 1     f(x): -0.500000
+    The algorithm reached approximately critical point after 1 iterations; the gradient norm (9.966301158124465e-7) is less than 0.001.
+    # 75    f(x): -0.500000 | ϵ: 0.00000562
+     | Initial f(x): -0.500000
+     | # 1     f(x): -0.500000
+    The algorithm reached approximately critical point after 1 iterations; the gradient norm (5.4875346930698466e-8) is less than 0.001.
+    # 100   f(x): -0.500000 | ϵ: 0.00000100
+    The value of the variable (ϵ) is smaller than or equal to its threshold (1.0e-6).
+    The algorithm performed a step with a change (6.5347623783315016e-9) less than 1.0e-6.
+
+where we now see that the subsolver always only requires one step. Note that since debug of an iteration is happening *after* a step, we see the sub solver run *before* the debug for an iteration number.

src/Manopt.jl

@@ -279,8 +279,6 @@ export get_proximal_map,
     set_gradient!,
     set_iterate!,
     set_manopt_parameter!,
-    set_manopt_parameter!,
-    set_manopt_parameter!,
     linearized_forward_operator,
     linearized_forward_operator!,
     adjoint_linearized_operator,
@@ -482,6 +480,7 @@ export DebugDualBaseChange, DebugDualBaseIterate, DebugDualChange, DebugDualIter
 export DebugDualResidual, DebugPrimalDualResidual, DebugPrimalResidual
 export DebugProximalParameter, DebugWarnIfCostIncreases
 export DebugGradient, DebugGradientNorm, DebugStepsize
+export DebugWhenActive
 export DebugWarnIfCostNotFinite, DebugWarnIfFieldNotFinite, DebugMessages
 #
 # Records - and access functions

src/plans/debug.jl

@@ -125,6 +125,16 @@ function show(io::IO, dg::DebugGroup)
     s = join(["$(di)" for di in dg.group], ", ")
     return print(io, "DebugGroup([$s])")
 end
+function set_manopt_parameter!(dg::DebugGroup, v::Val, args...)
+    for di in dg.group
+        set_manopt_parameter!(di, v, args...)
+    end
+    return dg
+end
+function set_manopt_parameter!(dg::DebugGroup, e::Symbol, args...)
+    set_manopt_parameter!(dg, Val(e), args...)
+    return dg
+end
 
 @doc raw"""
     DebugEvery <: DebugAction
@@ -154,6 +164,11 @@ function (d::DebugEvery)(p::AbstractManoptProblem, st::AbstractManoptSolverState
     elseif d.always_update
         d.debug(p, st, -1)
     end
+    # set activity for the next iterate in subsolvers
+    set_manopt_parameter!(
+        st, :SubState, :Debug, :active, !(i < 1) && (rem(i + 1, d.every) == 0)
+    )
+    return nothing
 end
 function show(io::IO, de::DebugEvery)
     return print(io, "DebugEvery($(de.debug), $(de.every), $(de.always_update))")
@@ -167,6 +182,15 @@ function status_summary(de::DebugEvery)
     end
     return "[$s, $(de.every)]"
 end
+function set_manopt_parameter!(de::DebugEvery, e::Symbol, args...)
+    set_manopt_parameter!(de, Val(e), args...)
+    return de
+end
+function set_manopt_parameter!(de::DebugEvery, args...)
+    set_manopt_parameter!(de.debug, args...)
+    return de
+end
+
 #
 # Special single ones
 #
@@ -256,13 +280,13 @@ print the current cost function value, see [`get_cost`](@ref).
 
 * `format` - (`"$prefix %f"`) format to print the output using sprintf and a prefix (see `long`).
 * `io` – (`stdout`) default steream to print the debug to.
-* `long` - (`false`) short form to set the format to `F(x):` (default) or `current cost: ` and the cost
+* `long` - (`false`) short form to set the format to `f(x):` (default) or `current cost: ` and the cost
 """
 mutable struct DebugCost <: DebugAction
     io::IO
     format::String
     function DebugCost(;
-        long::Bool=false, io::IO=stdout, format=long ? "current cost: %f" : "F(x): %f"
+        long::Bool=false, io::IO=stdout, format=long ? "current cost: %f" : "f(x): %f"
     )
         return new(io, format)
     end
@@ -576,6 +600,57 @@ end
 show(io::IO, ::DebugStoppingCriterion) = print(io, "DebugStoppingCriterion()")
 status_summary(::DebugStoppingCriterion) = ":Stop"
 
+@doc raw"""
+    DebugWhenActive <: DebugAction
+
+evaluate and print debug only if the active boolean is set.
+This can be set from outside and is for example triggered by [`DebugEvery`](@ref)
+on debugs on the subsolver.
+
+This method does not perform any print itself but relies on it's childrens print.
+
+For now, the main interaction is with [`DebugEvery`](@ref) which might activate or
+deactivate this debug
+
+# Fields
+* `always_update` – whether or not to call the order debugs with iteration `-1` in in active state
+* `active` – a boolean that can (de-)activated from outside to enable/disable debug
+
+# Constructor
+
+    DebugWhenActive(d::DebugAction, active=true, always_update=true)
+
+Initialise the DebugSubsolver.
+"""
+mutable struct DebugWhenActive <: DebugAction
+    debug::DebugAction
+    active::Bool
+    always_update::Bool
+    function DebugWhenActive(d::DebugAction, active::Bool=true, always_update::Bool=true)
+        return new(d, active, always_update)
+    end
+end
+function (dwa::DebugWhenActive)(p::AbstractManoptProblem, st::AbstractManoptSolverState, i)
+    if dwa.active
+        dwa.debug(p, st, i)
+    elseif dwa.always_update
+        dwa.debug(p, st, -1)
+    end
+end
+function show(io::IO, dwa::DebugWhenActive)
+    return print(io, "DebugWhenActive($(dwa.debug), $(dwa.active), $(dwa.always_update))")
+end
+function status_summary(dwa::DebugWhenActive)
+    return repr(dwa)
+end
+function set_manopt_parameter!(dwa::DebugWhenActive, v::Val, args...)
+    set_manopt_parameter!(dwa.debug, v, args...)
+    return dwa
+end
+function set_manopt_parameter!(dwa::DebugWhenActive, ::Val{:active}, v)
+    return dwa.active = v
+end
+
 @doc raw"""
     DebugTime()
 
@@ -808,6 +883,7 @@ given an array of `Symbol`s, `String`s [`DebugAction`](@ref)s and `Ints`
 
 * The symbol `:Stop` creates an entry of to display the stopping criterion at the end
   (`:Stop => DebugStoppingCriterion()`), for further symbols see [`DebugActionFactory`](@ref DebugActionFactory(::Symbol))
+* The symbol `:Subsolver` wraps all `dictionary` entries with [`DebugWhenActive`](@ref) that can be set from outside.
 * Tuples of a symbol and a string can be used to also specify a format, see [`DebugActionFactory`](@ref DebugActionFactory(::Tuple{Symbol,String}))
 * any string creates a [`DebugDivider`](@ref)
 * any [`DebugAction`](@ref) is directly included
@@ -836,7 +912,7 @@ It also adds the [`DebugStoppingCriterion`](@ref) to the `:Stop` entry of the di
 function DebugFactory(a::Array{<:Any,1})
     # filter out every
     group = Array{DebugAction,1}()
-    for s in filter(x -> !isa(x, Int) && x != :Stop, a) # filter ints and stop
+    for s in filter(x -> !isa(x, Int) && (x ∉ [:Stop, :Subsolver]), a) # filter ints and stop
         push!(group, DebugActionFactory(s))
     end
     dictionary = Dict{Symbol,DebugAction}()
@@ -849,8 +925,11 @@ function DebugFactory(a::Array{<:Any,1})
         end
         dictionary[:All] = debug
     end
-    if :Stop in a
-        dictionary[:Stop] = DebugStoppingCriterion()
+    (:Stop in a) && (dictionary[:Stop] = DebugStoppingCriterion())
+    if (:Subsolver in a)
+        for k in keys(dictionary)
+            dictionary[k] = DebugWhenActive(dictionary[k])
+        end
     end
     return dictionary
 end
@@ -887,8 +966,8 @@ Note that the Shortcut symbols should all start with a capital letter.
 * `:Time` creates a [`DebugTime`](@ref)
 * `:WarningMessages`creates a [`DebugMessages`](@ref)`(:Warning)`
 * `:InfoMessages`creates a [`DebugMessages`](@ref)`(:Info)`
-* `:ErrorMessages`creates a [`DebugMessages`](@ref)`(:Error)`
-* `:Messages`creates a [`DebugMessages`](@ref)`()` (i.e. the same as `:InfoMessages`)
+* `:ErrorMessages` creates a [`DebugMessages`](@ref)`(:Error)`
+* `:Messages` creates a [`DebugMessages`](@ref)`()` (i.e. the same as `:InfoMessages`)
 
 any other symbol creates a `DebugEntry(s)` to print the entry (o.:s) from the options.
 """

src/plans/difference_of_convex_plan.jl

@@ -129,7 +129,7 @@ function set_manopt_parameter!(ldc::LinearizedDCCost, ::Val{:p}, p)
     return ldc
 end
 function set_manopt_parameter!(ldc::LinearizedDCCost, ::Val{:X}, X)
-    ldc.Xk = X
+    ldc.Xk .= X
     return ldc
 end
 
@@ -203,7 +203,7 @@ function set_manopt_parameter!(ldcg::LinearizedDCGrad, ::Val{:p}, p)
     return ldcg
 end
 function set_manopt_parameter!(ldcg::LinearizedDCGrad, ::Val{:X}, X)
-    ldcg.Xk = X
+    ldcg.Xk .= X
     return ldcg
 end