diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 3822d29063..7e4221c8b6 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -11,7 +11,7 @@ jobs: runs-on: ${{ matrix.os }} strategy: matrix: - julia-version: ["1.6", "1.8", "~1.9.0-0"] + julia-version: ["1.6", "1.8", "1.9"] os: [ubuntu-latest, macOS-latest, windows-latest] steps: - uses: actions/checkout@v3 @@ -28,4 +28,4 @@ jobs: file: ./lcov.info name: codecov-umbrella fail_ci_if_error: false - if: ${{ matrix.julia-version == '1.8' && matrix.os =='ubuntu-latest' }} + if: ${{ matrix.os =='ubuntu-latest' }} diff --git a/Changelog.md b/Changelog.md index 3f22028c46..40de8af17c 100644 --- a/Changelog.md +++ b/Changelog.md @@ -5,6 +5,22 @@ All notable Changes to the Julia package `Manopt.jl` will be documented in this 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.4.21] - 22/05/2023 + +### Added + +* A `ManifoldCacheObjective` as a decorator for objectives to cache results of calls, + using LRU Caches as a weak dependency. For now this works with cost and gradient evaluations +* A `ManifoldCountObjective` as a decorator for objectives to enable counting of calls to for example the cost and the gradient +* adds a `return_objective` keyword, that switches the return of a solver to a tuple `(o, s)`, + where `o` is the (possibly decorated) objective, and `s` os the “classical” solver return (state or point). + This way the counted values can be accessed and the cache can be reused. +* change solvers on the mid level (form `solver(M, objective, p)`) to also accept decorated objectives + +### Changed +* Switch all Requires weak dependencies to actual weak dependencies starting in Julia 1.9 + + ## [0.4.20] - 11/05/2023 ### Changed @@ -33,7 +49,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added * A new interface of the form `alg(M, objective, p0)` to allow to reuse - objectives without creating `AbstractSolverStates` and calling `solve!`. This especially still allows for any decoration of the objective and/or the state using e.g. `debug=`, or `record=`. + objectives without creating `AbstractManoptSolverState`s and calling `solve!`. This especially still allows for any decoration of the objective and/or the state using e.g. `debug=`, or `record=`. ### Changed diff --git a/Project.toml b/Project.toml index 2a78a4d3d7..e2b7573f4e 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Manopt" uuid = "0fc0a36d-df90-57f3-8f93-d78a9fc72bb5" authors = ["Ronny Bergmann "] -version = "0.4.20" +version = "0.4.21" [deps] ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4" @@ -17,7 +17,6 @@ Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Requires = "ae029012-a4dd-5104-9daa-d747884805df" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" -StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" @@ -26,19 +25,32 @@ ColorSchemes = "3.5.0" ColorTypes = "0.9.1, 0.10, 0.11" Colors = "0.11.2, 0.12" DataStructures = "0.17, 0.18" +LRUCache = "1.4" ManifoldDiff = "0.2, 0.3" Manifolds = "0.8.57" ManifoldsBase = "0.14.4" Requires = "0.5, 1" -StaticArrays = "0.12, 1.0" julia = "1.6" +[extensions] +ManoptLRUCacheExt = "LRUCache" +ManoptLineSearchesExt = "LineSearches" +ManoptManifoldsExt = ["Manifolds"] +ManoptPlotsExt = "Plots" + [extras] ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" +LRUCache = "8ac3fa9e-de4c-5943-b1dc-09c6b5f20637" +LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255" Manifolds = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" -LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255" [targets] -test = ["Test", "ForwardDiff", "Manifolds", "Plots", "LineSearches"] +test = ["Test", "ForwardDiff", "Manifolds", "Plots", "LineSearches", "LRUCache"] + +[weakdeps] +LRUCache = "8ac3fa9e-de4c-5943-b1dc-09c6b5f20637" +LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255" +Manifolds = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" \ No newline at end of file diff --git a/Readme.md b/Readme.md index 26a92efaa5..1f1d3d7315 100644 --- a/Readme.md +++ b/Readme.md @@ -42,22 +42,6 @@ The following packages are related to `Manopt.jl` * [`Manifolds.jl`](https://juliamanifolds.github.io/Manifolds.jl/stable/) – a library of manifolds implemented using [`ManifoldsBase.jl`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/) :octocat: [GitHub repository](https://github.com/JuliaManifolds/Manifolds.jl) * [`ManifoldsDiff.jl`](https://juliamanifolds.github.io/ManifoldDiff.jl/stable/) – a package to use (Euclidean) AD tools on manifolds, that also provides several differentials and gradients. :octocat: [GitHub repository](https://github.com/JuliaManifolds/ManifoldDiff.jl) -## Further Packages & Links - -`Manopt.jl` belongs to the Manopt family: - -* [manopt.org](https://www.manopt.org) – The Matlab version of Manopt, see also their :octocat: [GitHub repository](https://github.com/NicolasBoumal/manopt) -* [pymanopt.org](https://www.pymanopt.org/) – The Python version of Manopt – providing also several AD backends, see also their :octocat: [GitHub repository](https://github.com/pymanopt/pymanopt) - -but there are also more packages providing tools on manifolds: - -* [Jax Geometry](https://bitbucket.org/stefansommer/jaxgeometry/src/main/) (Python/Jax) for differential geometry and stochastic dynamics with deep learning -* [Geomstats](https://geomstats.github.io) (Python with several backends) focusing on statistics and machine learning :octocat: [GitHub repository](https://github.com/geomstats/geomstats) -* [Geoopt](https://geoopt.readthedocs.io/en/latest/) (Python & PyTorch) – Riemannian ADAM & SGD. :octocat: [GitHub repository](https://github.com/geoopt/geoopt) -* [McTorch](https://github.com/mctorch/mctorch) (Python & PyToch) – Riemannian SGD, Adagrad, ASA & CG. -* [ROPTLIB](https://www.math.fsu.edu/~whuang2/papers/ROPTLIB.htm) (C++) a Riemannian OPTimization LIBrary :octocat: [GitHub repository](https://github.com/whuang08/ROPTLIB) -* [TF Riemopt](https://github.com/master/tensorflow-riemopt) (Python & TensorFlow) Riemannian optimization using TensorFlow - ## Citation If you use `Manopt.jl` in your work, please cite the following @@ -91,3 +75,21 @@ To refer to a certain version or the source code in general we recommend to cite for the most recent version or a corresponding version specific DOI, see [the list of all versions](https://zenodo.org/search?page=1&size=20&q=conceptrecid:%224290905%22&sort=-version&all_versions=True). Note that both citations are in [BibLaTeX](https://ctan.org/pkg/biblatex) format. + +## Further and Similar Packages & Links + +`Manopt.jl` belongs to the Manopt family: + +* [manopt.org](https://www.manopt.org) – The Matlab version of Manopt, see also their :octocat: [GitHub repository](https://github.com/NicolasBoumal/manopt) +* [pymanopt.org](https://www.pymanopt.org/) – The Python version of Manopt – providing also several AD backends, see also their :octocat: [GitHub repository](https://github.com/pymanopt/pymanopt) + +but there are also more packages providing tools on manifolds: + +* [Jax Geometry](https://bitbucket.org/stefansommer/jaxgeometry/src/main/) (Python/Jax) for differential geometry and stochastic dynamics with deep learning +* [Geomstats](https://geomstats.github.io) (Python with several backends) focusing on statistics and machine learning :octocat: [GitHub repository](https://github.com/geomstats/geomstats) +* [Geoopt](https://geoopt.readthedocs.io/en/latest/) (Python & PyTorch) – Riemannian ADAM & SGD. :octocat: [GitHub repository](https://github.com/geoopt/geoopt) +* [McTorch](https://github.com/mctorch/mctorch) (Python & PyToch) – Riemannian SGD, Adagrad, ASA & CG. +* [ROPTLIB](https://www.math.fsu.edu/~whuang2/papers/ROPTLIB.htm) (C++) a Riemannian OPTimization LIBrary :octocat: [GitHub repository](https://github.com/whuang08/ROPTLIB) +* [TF Riemopt](https://github.com/master/tensorflow-riemopt) (Python & TensorFlow) Riemannian optimization using TensorFlow + +Did you use `Manopt.jl` somewhere? Let us know! We'd love to collect those here as well. \ No newline at end of file diff --git a/_quarto.yml b/_quarto.yml index c7ae299f25..f52d3bd07d 100644 --- a/_quarto.yml +++ b/_quarto.yml @@ -23,4 +23,4 @@ format: variant: -raw_html wrap: none -jupyter: julia-1.8 +jupyter: julia-1.9 diff --git a/docs/Project.toml b/docs/Project.toml index 2413fd69ea..db6e72d8f3 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -7,8 +7,9 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000" Images = "916415d5-f1e6-5110-898d-aaa5f9f070e0" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" -LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +LRUCache = "8ac3fa9e-de4c-5943-b1dc-09c6b5f20637" LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" Manifolds = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e" ManifoldsBase = "3362f125-f0bb-47a3-aa74-596ffd7ef2fb" diff --git a/docs/make.jl b/docs/make.jl index 4afb7b97df..2258d9c784 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,5 +1,5 @@ using Documenter: DocMeta, HTML, MathJax3, deploydocs, makedocs -using Manopt, Manifolds, Plots +using LineSearches, LRUCache, Manopt, Manifolds, Plots generated_path = joinpath(@__DIR__, "src") base_url = "https://github.com/JuliaManifolds/Manopt.jl/blob/master/" @@ -31,6 +31,7 @@ makedocs(; "Get started: Optimize!" => "tutorials/Optimize!.md", "Speedup using Inplace computations" => "tutorials/InplaceGradient.md", "Use Automatic Differentiation" => "tutorials/AutomaticDifferentiation.md", + "Count and use a Cache" => "tutorials/CountAndCache.md", "Record values" => "tutorials/HowToRecord.md", "Do Contrained Optimization" => "tutorials/ConstrainedOptimization.md", "Do Geodesic Regression" => "tutorials/GeodesicRegression.md", diff --git a/docs/src/plans/objective.md b/docs/src/plans/objective.md index a6cc8e4ce5..2253391259 100644 --- a/docs/src/plans/objective.md +++ b/docs/src/plans/objective.md @@ -8,6 +8,7 @@ The Objective describes that actual cost function and all its properties. ```@docs AbstractManifoldObjective +AbstractDecoratedManifoldObjective ``` Which has two main different possibilities for its containing functions concerning the evaluation mode – not necessarily the cost, but for example gradient in an [`AbstractManifoldGradientObjective`](@ref). @@ -167,10 +168,34 @@ decorate_objective! Since single function calls, e.g. to the cost or the gradient, might be expensive, a simple cache objective exists as a decorator, that caches one cost value or gradient. -_This feature was just recently introduced in Manopt 0.4 and might still be a little unstable. -The `cache::Symbol=` keyword argument of the solvers might be extended or still change slightly for example._ +It can be activated/used with the `cache=` keyword argument available for every solver. ```@docs -SimpleCacheObjective -objective_cache_factory +Manopt.reset_counters! +Manopt.objective_cache_factory +``` + +#### A simple cache + +A first generic cache is always available, but it only caches one gradient and one cost function evaluation (for the same point). + +```@docs +SimpleManifoldCachedObjective +``` + +#### A Generic Cache + +For the more advanced cache, you need to implement some type of cache yourself, that provides a `get!` +and implement [`init_caches`](@ref). +This is for example provided if you load [`LRUCache.jl`](https://github.com/JuliaCollections/LRUCache.jl). Then you obtain + +```@docs +ManifoldCachedObjective +init_caches +``` + +### [Count Objective](@id ManifoldCountObjective) + +```@docs +ManifoldCountObjective ``` \ No newline at end of file diff --git a/docs/src/tutorials/CountAndCache.md b/docs/src/tutorials/CountAndCache.md new file mode 100644 index 0000000000..523475e19e --- /dev/null +++ b/docs/src/tutorials/CountAndCache.md @@ -0,0 +1,160 @@ +# How to Count and Cache Function Calls +Ronny Bergmann + +In this tutorial, we want to investigate the caching and counting (i.e. statistics) features +of [Manopt.jl](https://manoptjl.org). We will reuse the optimization tasks from the +introductionary tutorial [Get Started: Optimize!](https://manoptjl.org/stable/tutorials/Optimize!.html). + +## Introduction + +There are surely many ways to keep track for example of how often the cost function is called, +for example with a [functor](https://docs.julialang.org/en/v1/manual/methods/#Function-like-objects), as we used in an example in [How to Record Data](https://manoptjl.org/stable/tutorials/HowtoRecord.html) + +``` julia +mutable struct MyCost{I<:Integer} + count::I +end +MyCost() = MyCost{Int64}(0) +function (c::MyCost)(M, x) + c.count += 1 + # [ .. Actual implementation of the cost here ] +end +``` + +This still leaves a bit of work to the user, especially for tracking more than just the number of cost function evaluations. + +When the a function like objective or gradient is expensive to compute, it may make sense to cache its results. +Manopt.jl tries to minimize the number of repeated calls but sometimes they are necessary and harmless when the function is cheap to compute. +Caching of expensive function calls can for example be added using [Memoize.jl](https://github.com/JuliaCollections/Memoize.jl) by the user. +The approach in the solvers of [Manopt.jl](https://manoptjl.org) aims to simplify adding +both these capabilities on the level of calling a solver. + +## Technical Background + +The two ingdredients for a solver in [Manopt.jl](https://manoptjl.org) +are the [`AbstractManoptProblem`](@ref) and the [`AbstractManoptSolverState`](@ref), where the +former consists of the domain, that is the [manifold](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/types.html#The-AbstractManifold) and [`AbstractManifoldObjective`](@ref). + +Both recording and debug capabilities are implemented in a decorator pattern to the solver state. +They can be easily added using the `record=` and `debug=` in any solver call. +This pattern was recently extended, such that also the objective can be decorated. +This is how both caching and counting are implemented, as decorators of the [`AbstractManifoldObjective`](@ref) +and hence for example changing/extending the behaviour of a call to [`get_cost`](@ref). + +Let’s finish off the technical background by loading the necessary packages. +Besides [Manopt.jl](https://manoptjl.org) and [Manifolds.jl](https://juliamanifolds.github.io/Manifolds.jl/latest/) we also need +[LRUCaches.jl](https://github.com/JuliaCollections/LRUCache.jl) which are (since Julia 1.9) a weak dependency and provide +the *least recently used* strategy for our caches. + +``` julia +using Manopt, Manifolds, Random, LRUCache +``` + +## Counting + +We first define our task, the Riemannian Center of Mass from the [Get Started: Optimize!](https://manoptjl.org/stable/tutorials/Optimize!.html) tutorial. + +``` julia +n = 100 +σ = π / 8 +M = Sphere(2) +p = 1 / sqrt(2) * [1.0, 0.0, 1.0] +data = [exp(M, p, σ * rand(M; vector_at=p)) for i in 1:n]; +f(M, p) = sum(1 / (2 * n) * distance.(Ref(M), Ref(p), data) .^ 2) +grad_f(M, p) = sum(1 / n * grad_distance.(Ref(M), data, Ref(p))); +``` + +to now count how often the cost and the gradient are called, we use the `count=` keyword +argument that works in any solver to specify the elements of the objective whose calls we +want to count calls to. A full list is available in the documentation of the +[`AbstractManifoldObjective`](@ref). +To also see the result, we have to set `return_objective=true`. +This returns `(objective, p)` instead of just the solver result `p`. +We can further also set `return_state=true` to get even more information about the solver run. + +``` julia +gradient_descent(M, f, grad_f, data[1]; count=[:Cost, :Gradient], return_objective=true, return_state=true) +``` + + # Solver state for `Manopt.jl`s Gradient Descent + After 72 iterations + + ## Parameters + * retraction method: ExponentialRetraction() + + ## Stepsize + ArmijoLineseach() with keyword parameters + * initial_stepsize = 1.0 + * retraction_method = ExponentialRetraction() + * contraction_factor = 0.95 + * sufficient_decrease = 0.1 + + ## Stopping Criterion + Stop When _one_ of the following are fulfilled: + Max Iteration 200: not reached + |grad f| < 1.0e-9: reached + Overall: reached + This indicates convergence: Yes + + ## Statistics on function calls + * :Gradient : 217 + * :Cost : 298 + on a ManifoldGradientObjective{AllocatingEvaluation} + +And we see that statistics are shown in the end. To now also cache these calls, +we can use the `cache=` keyword argument. +Since now both the cache and the count “extend” the functionality of the objective, +the order is important: On the high-level interface, the `count` is treated first, which +means that only actual function calls and not cache look-ups are counted. +With the proper initialisation, you can use any caches here that support the +`get!(function, cache, key)!` update. All parts of the objective that can currently be cached are listed at [`ManifoldCachedObjective`](@ref). The solver call has a keyword `cache` that takes a tuple`(c, vs, n)` of three arguments, where `c` is a symbol for the type of cache, `vs` is a vector of symbols, which calls to cache and `n` is the size of the cache. If the last element is not provided, a suitable default (currently`n=10`) is used. + +Here we want to use `c=:LRU` caches for `vs=[Cost, :Gradient]` with a size of `n=25`. + +``` julia +r = gradient_descent(M, f, grad_f, data[1]; + count=[:Cost, :Gradient], + cache=(:LRU, [:Cost, :Gradient], 25), + return_objective=true, return_state=true) +``` + + # Solver state for `Manopt.jl`s Gradient Descent + After 72 iterations + + ## Parameters + * retraction method: ExponentialRetraction() + + ## Stepsize + ArmijoLineseach() with keyword parameters + * initial_stepsize = 1.0 + * retraction_method = ExponentialRetraction() + * contraction_factor = 0.95 + * sufficient_decrease = 0.1 + + ## Stopping Criterion + Stop When _one_ of the following are fulfilled: + Max Iteration 200: not reached + |grad f| < 1.0e-9: reached + Overall: reached + This indicates convergence: Yes + + ## Statistics on function calls + * :Gradient : 72 + * :Cost : 164 + on a ManifoldGradientObjective{AllocatingEvaluation} + +Since the default setup with [`ArmijoLinesearch`](@ref) needs the gradient and the +cost, and similarly the stopping criterion might (independently) evaluate the gradient, +the caching is quite helpful here. + +And of course also for this advanced return value of the solver, we can still access the +result as usual: + +``` julia +get_solver_result(r) +``` + + 3-element Vector{Float64}: + 0.7298774364923435 + 0.047665824852873 + 0.6819141418393224 diff --git a/ext/ManoptLRUCacheExt.jl b/ext/ManoptLRUCacheExt.jl new file mode 100644 index 0000000000..79659bd2b5 --- /dev/null +++ b/ext/ManoptLRUCacheExt.jl @@ -0,0 +1,93 @@ +module ManoptLRUCacheExt + +using Manopt +import Manopt: init_caches +using ManifoldsBase + +if isdefined(Base, :get_extension) + using LRUCache +else + # imports need to be relative for Requires.jl-based workflows: + # https://github.com/JuliaArrays/ArrayInterface.jl/pull/387 + using ..LRUCache +end + +# introduce LRU even as default. +function Manopt.init_caches( + M::AbstractManifold, caches::AbstractVector{<:Symbol}; kwargs... +) + return Manopt.init_caches(M, caches, LRU; kwargs...) +end + +""" + init_caches(caches, T::Type{LRU}; kwargs...) + +Given a vector of symbols `caches`, this function sets up the +`NamedTuple` of caches, where `T` is the type of cache to use. + +# Keyword arguments + +* `p` - (`rand(M)`) a point on a manifold, to both infere its type for keys and initialize caches +* `value` - (`0.0`) a value both typing and initialising number-caches, eg. for caching a cost. +* `X` - (`zero_vector(M, p)` a tangent vector at `p` to both type and initialize tangent vector caches +* `cache_size` - (`10`) a default cache size to use +* `cache_sizes` – (`Dict{Symbol,Int}()`) a dictionary of sizes for the `caches` to specify different (non-default) sizes +""" +function Manopt.init_caches( + M::AbstractManifold, + caches::AbstractVector{<:Symbol}, + ::Type{LRU}; + p::P=rand(M), + value::R=0.0, + X::T=zero_vector(M, p), + cache_size::Int=10, + cache_sizes::Dict{Symbol,Int}=Dict{Symbol,Int}(), +) where {P,R<:Real,T} + lru_caches = LRU[] + for c in caches + i = length(lru_caches) + m = get(cache_sizes, c, cache_size) + # Float cache, e.g. Cost + (c === :Cost) && push!(lru_caches, LRU{P,R}(; maxsize=m)) + # vectors – e.g. Constraints/EqCOnstraints/InEqCOnstraints + # (a) store whole vectors + (c === :EqualityConstraints) && push!(lru_caches, LRU{P,Vector{R}}(; maxsize=m)) + (c === :InequalityConstraints) && push!(lru_caches, LRU{P,Vector{R}}(; maxsize=m)) + (c === :Constraints) && push!(lru_caches, LRU{P,Vector{Vector{R}}}(; maxsize=m)) + # (b) store single entries, but with an point-index key + (c === :EqualityConstraint) && push!(lru_caches, LRU{Tuple{P,Int},R}(; maxsize=m)) + (c === :InequalityConstraint) && push!(lru_caches, LRU{Tuple{P,Int},R}(; maxsize=m)) + # Tangent Vector cache + # (a) the simple ones, like the gradient or the Hessian + (c === :Gradient) && push!(lru_caches, LRU{P,T}(; maxsize=m)) + (c === :SubGradient) && push!(lru_caches, LRU{P,T}(; maxsize=m)) + (c === :SubtrahendGradient) && push!(lru_caches, LRU{P,T}(; maxsize=m)) + # (b) indexed by point and vector + (c === :Hessian) && push!(lru_caches, LRU{Tuple{P,T},T}(; maxsize=m)) + (c === :Preconditioner) && push!(lru_caches, LRU{Tuple{P,T},T}(; maxsize=m)) + # (b) store tangent vectors of components, but with an point-index key + (c === :GradEqualityConstraint) && + push!(lru_caches, LRU{Tuple{P,Int},T}(; maxsize=m)) + (c === :GradInequalityConstraint) && + push!(lru_caches, LRU{Tuple{P,Int},T}(; maxsize=m)) + # For the (future) product tangent budle this might also be just Ts + (c === :GradEqualityConstraints) && push!(lru_caches, LRU{P,Vector{T}}(; maxsize=m)) + (c === :GradInequalityConstraints) && + push!(lru_caches, LRU{P,Vector{T}}(; maxsize=m)) + # (c === :StochasticGradient) + (c === :StochasticGradient) && push!(lru_caches, LRU{Tuple{P,Int},T}(; maxsize=m)) + (c === :StochasticGradients) && push!(lru_caches, LRU{P,Vector{T}}(; maxsize=m)) + # Point caches + # (b) proximal point - we have to again use (p, λ, i) as key + (c === :ProximalMap) && push!(lru_caches, LRU{Tuple{P,R,Int},P}(; maxsize=m)) + # None of the above matched -> unknown cache type + if length(lru_caches) == i #nothing pushed + error(""" + A cache for :$c seems to not be supported by LRU caches. + """) + end + end + return NamedTuple{Tuple(caches)}(lru_caches) +end + +end diff --git a/ext/ManoptLineSearchesExt.jl b/ext/ManoptLineSearchesExt.jl new file mode 100644 index 0000000000..807b63322e --- /dev/null +++ b/ext/ManoptLineSearchesExt.jl @@ -0,0 +1,70 @@ +module ManoptLineSearchesExt + +using Manopt +import Manopt: LineSearchesStepsize +using ManifoldsBase + +if isdefined(Base, :get_extension) + using LineSearches +else + # imports need to be relative for Requires.jl-based workflows: + # https://github.com/JuliaArrays/ArrayInterface.jl/pull/387 + using ..LineSearches +end + +function (cs::Manopt.LineSearchesStepsize)( + mp::AbstractManoptProblem, + s::AbstractManoptSolverState, + i::Int, + η=-get_gradient(s); + fp=get_cost(mp, get_iterate(s)), + kwargs..., +) + M = get_manifold(mp) + p = get_iterate(s) + X = get_gradient(s) + p_tmp = copy(M, p) + X_tmp = copy(M, p, X) + Y_tmp = copy(M, p, X) + f = get_cost_function(get_objective(mp)) + dphi_0 = real(inner(M, p, X, η)) + + # guess initial alpha + α0 = 1.0 + + # perform actual linesearch + + function ϕ(α) + retract!(M, p_tmp, p, η, α, cs.retraction_method) + return f(M, p_tmp) + end + function dϕ(α) + retract!(M, p_tmp, p, η, α, cs.retraction_method) + get_gradient!(mp, X_tmp, p_tmp) + vector_transport_to!(M, Y_tmp, p, η, p_tmp, cs.vector_transport_method) + return real(inner(M, p_tmp, Y_tmp, Y_tmp)) + end + function ϕdϕ(α) + # TODO: optimize? + retract!(M, p_tmp, p, η, α, cs.retraction_method) + get_gradient!(mp, X_tmp, p_tmp) + vector_transport_to!(M, Y_tmp, p, η, p_tmp, cs.vector_transport_method) + phi = f(M, p_tmp) + dphi = real(inner(M, p_tmp, Y_tmp, Y_tmp)) + return (phi, dphi) + end + + try + α, fp = cs.linesearch(ϕ, dϕ, ϕdϕ, α0, fp, dphi_0) + return α + catch ex + if isa(ex, LineSearches.LineSearchException) + # maybe indicate failure? + return zero(dphi_0) + else + rethrow(ex) + end + end +end + +end diff --git a/ext/ManoptManifoldsExt/ChambollePockManifolds.jl b/ext/ManoptManifoldsExt/ChambollePockManifolds.jl new file mode 100644 index 0000000000..3ff7ec25b9 --- /dev/null +++ b/ext/ManoptManifoldsExt/ChambollePockManifolds.jl @@ -0,0 +1,11 @@ +function Manopt.ChambollePockState( + M::AbstractManifold, + m::P, + n::Q, + p::P, + X::T; + N::AbstractManifold=TangentBundle(M), + kwargs..., +) where {P,Q,T} + return ChambollePockState(M, N, m, n, p, X; kwargs...) +end diff --git a/ext/ManoptManifoldsExt/ManoptManifoldsExt.jl b/ext/ManoptManifoldsExt/ManoptManifoldsExt.jl new file mode 100644 index 0000000000..10d4a816d3 --- /dev/null +++ b/ext/ManoptManifoldsExt/ManoptManifoldsExt.jl @@ -0,0 +1,33 @@ +module ManoptManifoldsExt + +using ManifoldsBase: exp, log, ParallelTransport, vector_transport_to +using Manopt +import Manopt: + max_stepsize, + prox_TV2, + grad_TV2, + alternating_gradient_descent, + alternating_gradient_descent!, + get_gradient, + get_gradient! +using LinearAlgebra: cholesky, det, diag, dot, Hermitian, qr, Symmetric, triu, I, Diagonal +import ManifoldsBase: copy, mid_point, mid_point! + +using ManifoldDiff: + adjoint_differential_shortest_geodesic_startpoint, + adjoint_differential_shortest_geodesic_endpoint + +if isdefined(Base, :get_extension) + using Manifolds +else + using ..Manifolds +end + +const NONMUTATINGMANIFOLDS = Union{Circle,PositiveNumbers,Euclidean{Tuple{}}} +include("manifold_functions.jl") +include("nonmutating_manifolds_functions.jl") +include("artificialDataFunctionsManifolds.jl") +include("ChambollePockManifolds.jl") +include("alternating_gradient.jl") + +end diff --git a/ext/ManoptManifoldsExt/alternating_gradient.jl b/ext/ManoptManifoldsExt/alternating_gradient.jl new file mode 100644 index 0000000000..ce7c208964 --- /dev/null +++ b/ext/ManoptManifoldsExt/alternating_gradient.jl @@ -0,0 +1,112 @@ + +@doc raw""" + X = get_gradient(M::ProductManifold, ago::ManifoldAlternatingGradientObjective, p) + get_gradient!(M::ProductManifold, P::ManifoldAlternatingGradientObjective, X, p) + +Evaluate all summands gradients at a point `p` on the `ProductManifold M` (in place of `X`) +""" +get_gradient(M::ProductManifold, ::ManifoldAlternatingGradientObjective, ::Any...) + +function get_gradient( + M::AbstractManifold, + mago::ManifoldAlternatingGradientObjective{AllocatingEvaluation,TC,<:AbstractVector}, + p, +) where {TC} + return ProductRepr([gi(M, p) for gi in mago.gradient!!]...) +end + +@doc raw""" + X = get_gradient(M::AbstractManifold, p::ManifoldAlternatingGradientObjective, p, k) + get_gradient!(M::AbstractManifold, p::ManifoldAlternatingGradientObjective, X, p, k) + +Evaluate one of the component gradients ``\operatorname{grad}f_k``, ``k∈\{1,…,n\}``, at `x` (in place of `Y`). +""" +function get_gradient( + M::ProductManifold, + mago::ManifoldAlternatingGradientObjective{AllocatingEvaluation,TC,<:Function}, + p, + k, +) where {TC} + return get_gradient(M, mago, p)[M, k] +end +function get_gradient!( + M::AbstractManifold, + X, + mago::ManifoldAlternatingGradientObjective{InplaceEvaluation,TC,<:AbstractVector}, + p, +) where {TC} + for (gi, Xi) in zip(mago.gradient!!, submanifold_components(M, X)) + gi(M, Xi, p) + end + return X +end + +function get_gradient!( + M::ProductManifold, + X, + mago::ManifoldAlternatingGradientObjective{AllocatingEvaluation,TC,<:Function}, + p, + k, +) where {TC} + copyto!(M[k], X, mago.gradient!!(M, p)[M, k]) + return X +end + +function alternating_gradient_descent( + M::ProductManifold, + f, + grad_f::Union{TgF,AbstractVector{<:TgF}}, + p=rand(M); + evaluation::AbstractEvaluationType=AllocatingEvaluation(), + kwargs..., +) where {TgF} + ago = ManifoldAlternatingGradientObjective(f, grad_f; evaluation=evaluation) + return alternating_gradient_descent(M, ago, p; evaluation=evaluation, kwargs...) +end +function alternating_gradient_descent( + M::ProductManifold, ago::ManifoldAlternatingGradientObjective, p; kwargs... +) + q = copy(M, p) + return alternating_gradient_descent!(M, ago, q; kwargs...) +end + +function alternating_gradient_descent!( + M::ProductManifold, + f, + grad_f::Union{TgF,AbstractVector{<:TgF}}, + p; + evaluation::AbstractEvaluationType=AllocatingEvaluation(), + kwargs..., +) where {TgF} + agmo = ManifoldAlternatingGradientObjective(f, grad_f; evaluation=evaluation) + return alternating_gradient_descent!(M, agmo, p; evaluation=evaluation, kwargs...) +end +function alternating_gradient_descent!( + M::ProductManifold, + agmo::ManifoldAlternatingGradientObjective, + p; + inner_iterations::Int=5, + stopping_criterion::StoppingCriterion=StopAfterIteration(100) | + StopWhenGradientNormLess(1e-9), + stepsize::Stepsize=default_stepsize(M, AlternatingGradientDescentState), + order_type::Symbol=:Linear, + order=collect(1:length(M.manifolds)), + retraction_method::AbstractRetractionMethod=default_retraction_method(M, typeof(p)), + kwargs..., +) + dagmo = Manopt.decorate_objective!(M, agmo; kwargs...) + dmp = DefaultManoptProblem(M, dagmo) + agds = AlternatingGradientDescentState( + M, + p; + inner_iterations=inner_iterations, + stopping_criterion=stopping_criterion, + stepsize=stepsize, + order_type=order_type, + order=order, + retraction_method=retraction_method, + ) + agds = Manopt.decorate_state!(agds; kwargs...) + Manopt.solve!(dmp, agds) + return Manopt.get_solver_return(get_objective(dmp), agds) +end diff --git a/ext/ManoptManifoldsExt/artificialDataFunctionsManifolds.jl b/ext/ManoptManifoldsExt/artificialDataFunctionsManifolds.jl new file mode 100644 index 0000000000..d1331eb42a --- /dev/null +++ b/ext/ManoptManifoldsExt/artificialDataFunctionsManifolds.jl @@ -0,0 +1,143 @@ + +function Manopt.artificial_S2_composite_bezier_curve() + M = Sphere(2) + d0 = [0.0, 0.0, 1.0] + d1 = [0.0, -1.0, 0.0] + d2 = [-1.0, 0.0, 0.0] + d3 = [0.0, 0.0, -1.0] + # + # control points - where b1- and b2- are constructed by the C1 condition + # + # We define three segments: 1 + b00 = d0 # also known as p0 + ξ0 = π / (8.0 * sqrt(2.0)) .* [1.0, -1.0, 0.0] # staring direction from d0 + b01 = exp(M, d0, ξ0) # b0+ + ξ1 = π / (4.0 * sqrt(2)) .* [1.0, 0.0, 1.0] + # b02 or b1- and b11 or b1+ are constructed by this vector with opposing sign + # to achieve a C1 curve + b02 = exp(M, d1, ξ1) + b03 = d1 + # 2 + b10 = d1 + b11 = exp(M, d1, -ξ1) # yields c1 condition + ξ2 = -π / (4 * sqrt(2)) .* [0.0, 1.0, -1.0] + b12 = exp(M, d2, ξ2) + b13 = d2 + # 3 + b20 = d2 + b21 = exp(M, d2, -ξ2) + ξ3 = π / (8.0 * sqrt(2)) .* [-1.0, 1.0, 0.0] + b22 = exp(M, d3, ξ3) + b23 = d3 + # hence the matrix of controlpoints for the curve reads + return [ + BezierSegment([b00, b01, b02, b03]), + BezierSegment([b10, b11, b12, b13]), + BezierSegment([b20, b21, b22, b23]), + ] +end + +function Manopt.artificial_S2_rotation_image( + pts::Int=64, rotations::Tuple{Float64,Float64}=(0.5, 0.5) +) + M = Sphere(2) + img = fill(zeros(3), pts, pts) + north = [1.0, 0.0, 0.0] + Rxy(a) = [cos(a) -sin(a) 0.0; sin(a) cos(a) 0.0; 0.0 0.0 1] + Rxz(a) = [cos(a) 0.0 -sin(a); 0.0 1.0 0.0; sin(a) 0.0 cos(a)] + for i in 1:pts + for j in 1:pts + x = i / pts * 2π * rotations[1] + y = j / pts * 2π * rotations[2] + img[i, j] = Rxy(x + y) * Rxz(x - y) * [0, 0, 1] + end + end + return img +end + +function Manopt.artificial_S2_lemniscate(p, t::Float64, a::Float64=π / 2.0) + M = Sphere(2) + tP = 2.0 * Float64(p[1] >= 0.0) - 1.0 # Take north or south pole + base = [0.0, 0.0, tP] + xc = a * (cos(t) / (sin(t)^2 + 1.0)) + yc = a * (cos(t) * sin(t) / (sin(t)^2 + 1.0)) + tV = vector_transport_to(M, base, [xc, yc, 0.0], p, ParallelTransport()) + return exp(M, p, tV) +end + +function Manopt.artificial_S2_whirl_image(pts::Int=64) + M = Sphere(2) + img = artificial_S2_rotation_image(pts, (0.5, 0.5)) + # Set WhirlPatches + sc = pts / 64 + patchSizes = floor.(sc .* [9, 9, 9, 9, 11, 11, 11, 15, 15, 15, 17, 21]) + patchCenters = + Integer.( + floor.( + sc .* + [[35, 7] [25, 41] [32, 25] [7, 60] [10, 5] [41, 58] [11, 41] [23, 56] [ + 38, 45 + ] [16, 28] [55, 42] [51, 16]], + ), + ) + patchSigns = [1, 1, -1, 1, -1, -1, 1, 1, -1, -1, 1, -1] + for i in 1:length(patchSizes) + pS = Integer(patchSizes[i]) + pSH = Integer(floor((patchSizes[i] - 1) / 2)) + pC = patchCenters[:, i] + pC = max.(Ref(1), pC .- pS) .+ pS + pSgn = patchSigns[i] + s = pS % 2 == 0 ? 1 : 0 + r = [pC[1] .+ ((-pSH):(pSH + s)), pC[2] .+ ((-pSH):(pSH + s))] + patch = artificial_S2_whirl_patch(pS) + if pSgn == -1 # opposite ? + patch = -patch + end + img[r...] = patch + end + return img +end + +function Manopt.artificial_SPD_image2(pts=64, fraction=0.66) + Zl = 4.0 * Matrix{Float64}(I, 3, 3) + # create a first matrix + α = 2.0 * π / 3 + β = π / 3 + B = [1.0 0.0 0.0; 0.0 cos(β) -sin(β); 0.0 sin(β) cos(β)] + A = [cos(α) -sin(α) 0.0; sin(α) cos(α) 0.0; 0.0 0.0 1.0] + Zo = Matrix(Symmetric(A * B * Diagonal([2.0, 4.0, 8.0]) * B' * A')) + # create a second matrix + α = -4.0 * π / 3 + β = -π / 3 + B = [1.0 0.0 0.0; 0.0 cos(β) -sin(β); 0.0 sin(β) cos(β)] + A = [cos(α) -sin(α) 0.0; sin(α) cos(α) 0.0; 0.0 0.0 1.0] + Zt = A * B * Diagonal([8.0 / sqrt(2.0), 8.0, sqrt(2.0)]) * B' * A' + data = fill(Matrix{Float64}(I, 3, 3), pts, pts) + M = SymmetricPositiveDefinite(3) + for row in 1:pts + for col in 1:pts + # (a) from Zo a part to Zt + C = Zo + if (row > 1) # in X direction + C = exp( + M, + C, + log(M, C, Zt), + (row - 1) / (2 * (pts - 1)) + ((row > fraction * pts) ? 1 / 2 : 0.0), + ) + end + if (col > 1) # and then in Y direction + C = exp( + M, + C, + vector_transport_to( + M, Symmetric(Zo), log(M, Zo, Zl), Symmetric(C), ParallelTransport() + ), + (col - 1.0) / (pts - 1), + ) + end + data[row, col] = C + end + end + return data +end diff --git a/ext/ManoptManifoldsExt/manifold_functions.jl b/ext/ManoptManifoldsExt/manifold_functions.jl new file mode 100644 index 0000000000..d02fa3d71a --- /dev/null +++ b/ext/ManoptManifoldsExt/manifold_functions.jl @@ -0,0 +1,68 @@ + +""" + max_stepsize(M::TangentBundle, p) + +Tangent bundle has injectivity radius of either infinity (for flat manifolds) or 0 +(for non-flat manifolds). This makes a guess of what a reasonable maximum stepsize +on a tangent bundle might be. +""" +function max_stepsize(M::TangentBundle, p) + return max_stepsize(M.manifold, p[M, :point]) +end + +""" + mid_point(M, p, q, x) + mid_point!(M, y, p, q, x) + +Compute the mid point between `p` and `q`. If there is more than one mid point +of (not neccessarily minimizing) geodesics (e.g. on the sphere), the one nearest +to `x` is returned (in place of `y`). +""" +mid_point(M::AbstractManifold, p, q, ::Any) = mid_point(M, p, q) +mid_point!(M::AbstractManifold, y, p, q, ::Any) = mid_point!(M, y, p, q) + +function mid_point(M::Circle, p, q, x) + if distance(M, p, q) ≈ π + X = 0.5 * log(M, p, q) + Y = log(M, p, x) + return exp(M, p, (sign(X) == sign(Y) ? 1 : -1) * X) + end + return mid_point(M, p, q) +end + +function mid_point(M::Sphere, p, q, x) + if isapprox(M, p, -q) + X = log(M, p, x) / distance(M, p, x) * π + else + X = log(M, p, q) + end + return exp(M, p, 0.5 * X) +end +function mid_point!(M::Sphere, y, p, q, x) + if isapprox(M, p, -q) + X = log(M, p, x) / distance(M, p, x) * π + else + X = log(M, p, q) + end + y .= exp(M, p, 0.5 * X) + return y +end + +function prox_TV2(::Euclidean, λ, pointTuple::Tuple{T,T,T}, p::Int=1) where {T} + w = [1.0, -2.0, 1.0] + x = [pointTuple...] + if p == 1 # Example 3.2 in Bergmann, Laus, Steidl, Weinmann, 2014. + m = min.(Ref(λ), abs.(x .* w) / (dot(w, w))) + s = sign.(sum(x .* w)) + return x .- m .* s .* w + elseif p == 2 # Theorem 3.6 ibd. + t = λ * sum(x .* w) / (1 + λ * dot(w, w)) + return x .- t .* w + else + throw( + ErrorException( + "Proximal Map of TV2(Euclidean,λ,pT,p) not implemented for p=$(p) (requires p=1 or 2)", + ), + ) + end +end diff --git a/src/functions/nonmutating_manifolds_functions.jl b/ext/ManoptManifoldsExt/nonmutating_manifolds_functions.jl similarity index 97% rename from src/functions/nonmutating_manifolds_functions.jl rename to ext/ManoptManifoldsExt/nonmutating_manifolds_functions.jl index 1b68d4aa55..5c771bd51b 100644 --- a/src/functions/nonmutating_manifolds_functions.jl +++ b/ext/ManoptManifoldsExt/nonmutating_manifolds_functions.jl @@ -25,8 +25,8 @@ function grad_TV2(M::NONMUTATINGMANIFOLDS, q, p::Int=1) return X end function prox_TV2(::NONMUTATINGMANIFOLDS, λ, pointTuple::Tuple{T,T,T}, p::Int=1) where {T} - w = @SVector [1.0, -2.0, 1.0] - x = SVector(pointTuple) + w = [1.0, -2.0, 1.0] + x = [pointTuple...] if p == 1 # Theorem 3.5 in Bergmann, Laus, Steidl, Weinmann, 2014. sr_dot_xw = sym_rem(sum(x .* w)) m = min(λ, abs(sr_dot_xw) / (dot(w, w))) diff --git a/ext/ManoptPlotsExt/ManoptPlotsExt.jl b/ext/ManoptPlotsExt/ManoptPlotsExt.jl new file mode 100644 index 0000000000..874cc88dee --- /dev/null +++ b/ext/ManoptPlotsExt/ManoptPlotsExt.jl @@ -0,0 +1,14 @@ +module ManoptPlotsExt + +using Manopt +using Printf + +if isdefined(Base, :get_extension) + using Plots +else + # imports need to be relative for Requires.jl-based workflows: + # https://github.com/JuliaArrays/ArrayInterface.jl/pull/387 + using ..Plots +end +include("check_plots.jl") +end diff --git a/src/helpers/check_plots.jl b/ext/ManoptPlotsExt/check_plots.jl similarity index 58% rename from src/helpers/check_plots.jl rename to ext/ManoptPlotsExt/check_plots.jl index 5c70880ce7..504ff419ef 100644 --- a/src/helpers/check_plots.jl +++ b/ext/ManoptPlotsExt/check_plots.jl @@ -1,14 +1,5 @@ -""" -plot_slope(x, y; slope=2, line_base=0, a=0, b=2.0, i=1,j=length(x)) -plot the result from the error check functions, e.g. -[`check_gradient`](@ref), [`check_differential`](@ref), [`check_Hessian`](@ref) -on data `x,y` with two comparison lines - -1) `line_base` + t`slope` as the global slope the plot should have -2) `a` + `b*t` on the interval [`x[i]`, `x[j]`] for some (best fitting) comparison slope -""" -function plot_slope(x, y; slope=2, line_base=0, a=0, b=2.0, i=1, j=length(x)) +function Manopt.plot_slope(x, y; slope=2, line_base=0, a=0, b=2.0, i=1, j=length(x)) fig = plot( x, y; diff --git a/src/Manopt.jl b/src/Manopt.jl index 3e9b232fdc..85b6316b33 100644 --- a/src/Manopt.jl +++ b/src/Manopt.jl @@ -1,22 +1,21 @@ @doc raw""" 🏔️ Manopt.jl – Optimization on Manifolds in Julia. -- 📚 Documentation: [manoptjl.org](https://manoptjl.org) -- 📦 Repository: [github.com/JuliaManifolds/Manopt.jl](https://github.com/JuliaManifolds/Manopt.jl) -- 💬 Discussions: [github.com/JuliaManifolds/Manopt.jl/discussions](https://github.com/JuliaManifolds/Manopt.jl/discussions) -- 🎯 Issues: [github.com/JuliaManifolds/Manopt.jl/issues](https://github.com/JuliaManifolds/Manopt.jl/issues) +* 📚 Documentation: [manoptjl.org](https://manoptjl.org) +* 📦 Repository: [github.com/JuliaManifolds/Manopt.jl](https://github.com/JuliaManifolds/Manopt.jl) +* 💬 Discussions: [github.com/JuliaManifolds/Manopt.jl/discussions](https://github.com/JuliaManifolds/Manopt.jl/discussions) +* 🎯 Issues: [github.com/JuliaManifolds/Manopt.jl/issues](https://github.com/JuliaManifolds/Manopt.jl/issues) """ module Manopt import Base: &, copy, getindex, identity, setindex!, show, | import LinearAlgebra: reflect! -import ManifoldsBase: mid_point, mid_point! using ColorSchemes using ColorTypes using Colors using DataStructures: CircularBuffer, capacity, length, push!, size using Dates: Millisecond, Nanosecond, Period, canonicalize, value -using LinearAlgebra: Diagonal, I, eigen, eigvals, tril +using LinearAlgebra: Diagonal, I, eigen, eigvals, tril, Symmetric, dot, cholesky using ManifoldDiff: adjoint_Jacobi_field, adjoint_Jacobi_field!, @@ -99,6 +98,8 @@ using ManifoldsBase: log, log!, manifold_dimension, + mid_point, + mid_point!, norm, number_eltype, power_dimensions, @@ -120,10 +121,9 @@ using ManifoldsBase: ℝ using Markdown using Printf -using Random: shuffle! +using Random: shuffle!, rand, randperm using Requires using SparseArrays -using StaticArrays using Statistics: cor, cov, mean, std include("plans/plan.jl") @@ -134,9 +134,11 @@ include("functions/costs.jl") include("functions/differentials.jl") include("functions/gradients.jl") include("functions/proximal_maps.jl") +include("functions/manifold_functions.jl") # solvers general framework include("solvers/solver.jl") # specific solvers +include("solvers/alternating_gradient_descent.jl") include("solvers/augmented_Lagrangian_method.jl") include("solvers/ChambollePock.jl") include("solvers/conjugate_gradient_descent.jl") @@ -161,70 +163,46 @@ include("solvers/record_solver.jl") include("helpers/checks.jl") include("helpers/errorMeasures.jl") include("helpers/exports/Asymptote.jl") +include("helpers/LineSearchesTypes.jl") include("data/artificialDataFunctions.jl") include("deprecated.jl") function __init__() - @require Manifolds = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e" begin - using .Manifolds: - Circle, - Euclidean, - Grassmann, - GroupManifold, - Hyperbolic, - PositiveNumbers, - ProductManifold, - Rotations, - SymmetricPositiveDefinite, - Stiefel, - Sphere, - TangentBundle, - TangentSpaceAtPoint, - FixedRankMatrices, - SVDMPoint, - UMVTVector, - ArrayPowerRepresentation, - ProductRepr, - submanifold_components, - sym_rem, - mean - import Random: rand, randperm - import ManifoldsBase: copy - using LinearAlgebra: cholesky, det, diag, dot, Hermitian, qr, Symmetric, triu - # adaptions for Nonmutating manifolds - const NONMUTATINGMANIFOLDS = Union{Circle,PositiveNumbers,Euclidean{Tuple{}}} - include("functions/manifold_functions.jl") - include("functions/nonmutating_manifolds_functions.jl") - include("plans/alternating_gradient_plan.jl") - include("solvers/alternating_gradient_descent.jl") - export mid_point, mid_point!, reflect, reflect! - export AlternatingGradientDescentState - export AlternatingGradient - export alternating_gradient_descent, alternating_gradient_descent! - end - @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin - using .Plots - include("helpers/check_plots.jl") - end - @require LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255" begin - using .LineSearches - include("ext/LineSearchesExt.jl") + # + # Requires fallback for Julia < 1.9 + # + @static if !isdefined(Base, :get_extension) + @require Manifolds = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e" begin + include("../ext/ManoptManifoldsExt/ManoptManifoldsExt.jl") + end + @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin + include("../ext/ManoptPlotsExt/ManoptPlotsExt.jl") + end + @require LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255" begin + include("../ext/ManoptLineSearchesExt.jl") + end + @require LRUCache = "8ac3fa9e-de4c-5943-b1dc-09c6b5f20637" begin + include("../ext/ManoptLRUCacheExt.jl") + end end return nothing end # # General export ℝ, ℂ, &, | +export mid_point, mid_point!, reflect, reflect! # # Problems export AbstractManoptProblem, DefaultManoptProblem, TwoManifoldProblem # # Objectives -export AbstractManifoldGradientObjective, +export AbstractDecoratedManifoldObjective, + AbstractManifoldGradientObjective, AbstractManifoldCostObjective, AbstractManifoldObjective, AbstractPrimalDualManifoldObjective, ConstrainedManifoldObjective, + ManifoldCountObjective, NonlinearLeastSquaresObjective, ManifoldAlternatingGradientObjective, ManifoldCostGradientObjective, @@ -238,7 +216,8 @@ export AbstractManifoldGradientObjective, ManifoldSubgradientObjective, PrimalDualManifoldObjective, PrimalDualManifoldSemismoothNewtonObjective, - SimpleCacheObjective + SimpleManifoldCachedObjective, + ManifoldCachedObjective # # Evaluation & Problems - old export AbstractEvaluationType, AllocatingEvaluation, InplaceEvaluation, evaluation_type @@ -248,6 +227,7 @@ export AbstractGradientSolverState, AbstractHessianSolverState, AbstractManoptSolverState, AbstractPrimalDualSolverState, + AlternatingGradientDescentState, AugmentedLagrangianMethodState, ChambollePockState, ConjugateGradientDescentState, @@ -270,6 +250,7 @@ export AbstractGradientSolverState, export FrankWolfeCost, FrankWolfeGradient export NelderMeadSimplex +export AlternatingGradient # # Accessors and helpers for AbstractManoptSolverState export default_stepsize @@ -356,7 +337,9 @@ export DirectionUpdateRule, ConjugateGradientBealeRestart # # Solvers -export augmented_Lagrangian_method, +export alternating_gradient_descent, + alternating_gradient_descent!, + augmented_Lagrangian_method, augmented_Lagrangian_method!, ChambollePock, ChambollePock!, @@ -395,7 +378,7 @@ export augmented_Lagrangian_method, trust_regions! # Solver helpers export decorate_state!, decorate_objective! -export initialize_solver!, step_solver!, get_solver_result, get_solver_return, stop_solver! +export initialize_solver!, step_solver!, get_solver_result, stop_solver! export solve! export ApproxHessianFiniteDifference, ApproxHessianSymmetricRankOne, ApproxHessianBFGS export update_hessian!, update_hessian_basis! @@ -512,6 +495,9 @@ export RecordPrimalBaseChange, export RecordDualBaseChange, RecordDualBaseIterate, RecordDualChange, RecordDualIterate export RecordProximalParameter # +# Count +export get_count, reset_counters! +# # Helpers export check_gradient, check_differential, check_Hessian end diff --git a/src/data/artificialDataFunctions.jl b/src/data/artificialDataFunctions.jl index 377b4b13dc..4db50ed9a6 100644 --- a/src/data/artificialDataFunctions.jl +++ b/src/data/artificialDataFunctions.jl @@ -164,9 +164,13 @@ function artificial_S1_signal(x::Real) end return y end + +function artificial_S2_whirl_image end + @doc raw""" - artificial_S2_whirl_image([pts=64]) -generate an artificial image of data on the 2 sphere, + artificial_S2_whirl_image([pts::Int=64]) + +Generate an artificial image of data on the 2 sphere, # Arguments * `pts` – (`64`) size of the image in `pts`$\times$`pts` pixel. @@ -181,41 +185,14 @@ This example dataset was used in the numerical example in Section 5.5 of It is based on [`artificial_S2_rotation_image`](@ref) extended by small whirl patches. """ -function artificial_S2_whirl_image(pts::Int=64) - M = Sphere(2) - img = artificial_S2_rotation_image(pts, (0.5, 0.5)) - # Set WhirlPatches - sc = pts / 64 - patchSizes = floor.(sc .* [9, 9, 9, 9, 11, 11, 11, 15, 15, 15, 17, 21]) - patchCenters = - Integer.( - floor.( - sc .* - [[35, 7] [25, 41] [32, 25] [7, 60] [10, 5] [41, 58] [11, 41] [23, 56] [ - 38, 45 - ] [16, 28] [55, 42] [51, 16]], - ), - ) - patchSigns = [1, 1, -1, 1, -1, -1, 1, 1, -1, -1, 1, -1] - for i in 1:length(patchSizes) - pS = Integer(patchSizes[i]) - pSH = Integer(floor((patchSizes[i] - 1) / 2)) - pC = patchCenters[:, i] - pC = max.(Ref(1), pC .- pS) .+ pS - pSgn = patchSigns[i] - s = pS % 2 == 0 ? 1 : 0 - r = [pC[1] .+ ((-pSH):(pSH + s)), pC[2] .+ ((-pSH):(pSH + s))] - patch = artificial_S2_whirl_patch(pS) - if pSgn == -1 # opposite ? - patch = -patch - end - img[r...] = patch - end - return img -end +artificial_S2_whirl_image(::Int) + +function artificial_S2_rotation_image end + @doc raw""" artificial_S2_rotation_image([pts=64, rotations=(.5,.5)]) -creates an image with a rotation on each axis as a parametrization. + +Create an image with a rotation on each axis as a parametrization. # Optional Parameters * `pts` – (`64`) number of pixels along one dimension @@ -230,23 +207,7 @@ This dataset was used in the numerical example of Section 5.1 of > arxiv: [1506.02409](https://arxiv.org/abs/1506.02409) """ -function artificial_S2_rotation_image( - pts::Int=64, rotations::Tuple{Float64,Float64}=(0.5, 0.5) -) - M = Sphere(2) - img = fill(zeros(3), pts, pts) - north = [1.0, 0.0, 0.0] - Rxy(a) = [cos(a) -sin(a) 0.0; sin(a) cos(a) 0.0; 0.0 0.0 1] - Rxz(a) = [cos(a) 0.0 -sin(a); 0.0 1.0 0.0; sin(a) 0.0 cos(a)] - for i in 1:pts - for j in 1:pts - x = i / pts * 2π * rotations[1] - y = j / pts * 2π * rotations[2] - img[i, j] = Rxy(x + y) * Rxz(x - y) * [0, 0, 1] - end - end - return img -end +artificial_S2_rotation_image() @doc raw""" artificial_S2_whirl_patch([pts=5]) @@ -276,6 +237,9 @@ function artificial_S2_whirl_patch(pts::Int=5) end return patch end + +function artificial_S2_composite_bezier_curve end + @doc raw""" artificial_S2_composite_bezier_curve() @@ -314,43 +278,7 @@ This example was used within minimization of acceleration of the paper > arxiv: [1807.10090](https://arxiv.org/abs/1807.10090) """ -function artificial_S2_composite_bezier_curve() - M = Sphere(2) - d0 = [0.0, 0.0, 1.0] - d1 = [0.0, -1.0, 0.0] - d2 = [-1.0, 0.0, 0.0] - d3 = [0.0, 0.0, -1.0] - # - # control points - where b1- and b2- are constructed by the C1 condition - # - # We define three segments: 1 - b00 = d0 # also known as p0 - ξ0 = π / (8.0 * sqrt(2.0)) .* [1.0, -1.0, 0.0] # staring direction from d0 - b01 = exp(M, d0, ξ0) # b0+ - ξ1 = π / (4.0 * sqrt(2)) .* [1.0, 0.0, 1.0] - # b02 or b1- and b11 or b1+ are constructed by this vector with opposing sign - # to achieve a C1 curve - b02 = exp(M, d1, ξ1) - b03 = d1 - # 2 - b10 = d1 - b11 = exp(M, d1, -ξ1) # yields c1 condition - ξ2 = -π / (4 * sqrt(2)) .* [0.0, 1.0, -1.0] - b12 = exp(M, d2, ξ2) - b13 = d2 - # 3 - b20 = d2 - b21 = exp(M, d2, -ξ2) - ξ3 = π / (8.0 * sqrt(2)) .* [-1.0, 1.0, 0.0] - b22 = exp(M, d3, ξ3) - b23 = d3 - # hence the matrix of controlpoints for the curve reads - return [ - BezierSegment([b00, b01, b02, b03]), - BezierSegment([b10, b11, b12, b13]), - BezierSegment([b20, b21, b22, b23]), - ] -end +artificial_S2_composite_bezier_curve() @doc raw""" artificial_SPD_image([pts=64, stepsize=1.5]) @@ -392,6 +320,9 @@ function artificial_SPD_image(pts::Int=64, stepsize=1.5) end return data end + +function artificial_SPD_image2 end + @doc raw""" artificial_SPD_image2([pts=64, fraction=.66]) @@ -407,65 +338,23 @@ This data set was introduced in the numerical examples of Section of > arxiv: [1512.02814](https://arxiv.org/abs/1512.02814) """ -function artificial_SPD_image2(pts=64, fraction=0.66) - Zl = 4.0 * Matrix{Float64}(I, 3, 3) - # create a first matrix - α = 2.0 * π / 3 - β = π / 3 - B = [1.0 0.0 0.0; 0.0 cos(β) -sin(β); 0.0 sin(β) cos(β)] - A = [cos(α) -sin(α) 0.0; sin(α) cos(α) 0.0; 0.0 0.0 1.0] - Zo = Matrix(Symmetric(A * B * Diagonal([2.0, 4.0, 8.0]) * B' * A')) - # create a second matrix - α = -4.0 * π / 3 - β = -π / 3 - B = [1.0 0.0 0.0; 0.0 cos(β) -sin(β); 0.0 sin(β) cos(β)] - A = [cos(α) -sin(α) 0.0; sin(α) cos(α) 0.0; 0.0 0.0 1.0] - Zt = A * B * Diagonal([8.0 / sqrt(2.0), 8.0, sqrt(2.0)]) * B' * A' - data = fill(Matrix{Float64}(I, 3, 3), pts, pts) - M = SymmetricPositiveDefinite(3) - for row in 1:pts - for col in 1:pts - # (a) from Zo a part to Zt - C = Zo - if (row > 1) # in X direction - C = exp( - M, - C, - log(M, C, Zt), - (row - 1) / (2 * (pts - 1)) + ((row > fraction * pts) ? 1 / 2 : 0.0), - ) - end - if (col > 1) # and then in Y direction - C = exp( - M, - C, - vector_transport_to( - M, Symmetric(Zo), log(M, Zo, Zl), Symmetric(C), ParallelTransport() - ), - (col - 1.0) / (pts - 1), - ) - end - data[row, col] = C - end - end - return data -end +artificial_SPD_image2(pts, fraction) @doc raw""" - artificial_S2_lemniscate(p [,pts=128,a=π/2,interval=[0,2π]) + artificial_S2_lemniscate(p, t::Float64; a::Float64=π/2) -generate a Signal on the [Sphere](https://juliamanifolds.github.io/Manifolds.jl/stable/manifolds/sphere.html) $\mathbb S^2$ by creating the -[Lemniscate of Bernoulli](https://en.wikipedia.org/wiki/Lemniscate_of_Bernoulli) -in the tangent space of `p` sampled at `pts` points and use `exp` to get a -signal on the [Sphere](https://juliamanifolds.github.io/Manifolds.jl/stable/manifolds/sphere.html). +Generate a point from the signal on the [`Sphere`](https://juliamanifolds.github.io/Manifolds.jl/stable/manifolds/sphere.html) $\mathbb S^2$ by +creating the [Lemniscate of Bernoulli](https://en.wikipedia.org/wiki/Lemniscate_of_Bernoulli) +in the tangent space of `p` sampled at `t` and use èxp` to obtain a point on +the [`Sphere`](https://juliamanifolds.github.io/Manifolds.jl/stable/manifolds/sphere.html). # Input * `p` – the tangent space the Lemniscate is created in -* `pts` – (`128`) number of points to sample the Lemniscate -* `a` – (`π/2`) defines a half axis of the Lemniscate to cover a +* `t` – value to sample the Lemniscate at + +# Optional Values + * `a` – (`π/2`) defines a half axis of the Lemniscate to cover a half sphere. -* `interval` – (`[0,2*π]`) range to sample the lemniscate at, the default value - refers to one closed curve This dataset was used in the numerical example of Section 5.1 of @@ -476,26 +365,23 @@ This dataset was used in the numerical example of Section 5.1 of > arxiv: [1506.02409](https://arxiv.org/abs/1506.02409) """ -function artificial_S2_lemniscate( - p, pts::Integer=128, a::Float64=π / 2.0, interval::Array{Float64,1}=[0.0, 2.0 * π] -) - return artificial_S2_lemniscate.(Ref(p), range(interval[1], interval[2]; length=pts), a) -end +artificial_S2_lemniscate(p, t::Float64, a::Float64=π / 2.0) + @doc raw""" - artificial_S2_lemniscate(p,t; a=π/2) + artificial_S2_lemniscate(p [,pts=128,a=π/2,interval=[0,2π]) -generate a point from the signal on the [Sphere](https://juliamanifolds.github.io/Manifolds.jl/stable/manifolds/sphere.html) $\mathbb S^2$ by -creating the [Lemniscate of Bernoulli](https://en.wikipedia.org/wiki/Lemniscate_of_Bernoulli) -in the tangent space of `p` sampled at `t` and use èxp` to obtain a point on -the [Sphere](https://juliamanifolds.github.io/Manifolds.jl/stable/manifolds/sphere.html). +Generate a Signal on the [`Sphere`](https://juliamanifolds.github.io/Manifolds.jl/stable/manifolds/sphere.html) $\mathbb S^2$ by creating the +[Lemniscate of Bernoulli](https://en.wikipedia.org/wiki/Lemniscate_of_Bernoulli) +in the tangent space of `p` sampled at `pts` points and use `exp` to get a +signal on the [`Sphere`](https://juliamanifolds.github.io/Manifolds.jl/stable/manifolds/sphere.html). # Input * `p` – the tangent space the Lemniscate is created in -* `t` – value to sample the Lemniscate at - -# Optional Values - * `a` – (`π/2`) defines a half axis of the Lemniscate to cover a +* `pts` – (`128`) number of points to sample the Lemniscate +* `a` – (`π/2`) defines a half axis of the Lemniscate to cover a half sphere. +* `interval` – (`[0,2*π]`) range to sample the lemniscate at, the default value + refers to one closed curve This dataset was used in the numerical example of Section 5.1 of @@ -505,13 +391,9 @@ This dataset was used in the numerical example of Section 5.1 of > doi: [10.1137/15M101988X](https://dx.doi.org/10.1137/15M101988X) > arxiv: [1506.02409](https://arxiv.org/abs/1506.02409) - """ -function artificial_S2_lemniscate(p, t::Float64, a::Float64=π / 2.0) - M = Sphere(2) - tP = 2.0 * Float64(p[1] >= 0.0) - 1.0 # Take north or south pole - base = [0.0, 0.0, tP] - xc = a * (cos(t) / (sin(t)^2 + 1.0)) - yc = a * (cos(t) * sin(t) / (sin(t)^2 + 1.0)) - tV = vector_transport_to(M, base, [xc, yc, 0.0], p, ParallelTransport()) - return exp(M, p, tV) +""" +function artificial_S2_lemniscate( + p, pts::Integer=128, a::Float64=π / 2.0, interval::Array{Float64,1}=[0.0, 2.0 * π] +) + return artificial_S2_lemniscate.(Ref(p), range(interval[1], interval[2]; length=pts), a) end diff --git a/src/deprecated.jl b/src/deprecated.jl index 1d03d80b1b..53d39698b1 100644 --- a/src/deprecated.jl +++ b/src/deprecated.jl @@ -1,2 +1,4 @@ Base.@deprecate_binding HeestenesStiefelCoefficient HestenesStiefelCoefficient export HeestenesStiefelCoefficient +Base.@deprecate_binding SimpleCacheObjective SimpleManifoldCachedObjective +export SimpleCacheObjective diff --git a/src/functions/manifold_functions.jl b/src/functions/manifold_functions.jl index 8d3c9b0dec..e31d7dcc8f 100644 --- a/src/functions/manifold_functions.jl +++ b/src/functions/manifold_functions.jl @@ -1,71 +1,3 @@ - -""" - max_stepsize(M::TangentBundle, p) - -Tangent bundle has injectivity radius of either infinity (for flat manifolds) or 0 -(for non-flat manifolds). This makes a guess of what a reasonable maximum stepsize -on a tangent bundle might be. -""" -function max_stepsize(M::TangentBundle, p) - return max_stepsize(M.manifold, p[M, :point]) -end - -""" - mid_point(M, p, q, x) - mid_point!(M, y, p, q, x) - -Compute the mid point between `p` and `q`. If there is more than one mid point -of (not neccessarily minimizing) geodesics (e.g. on the sphere), the one nearest -to `x` is returned (in place of `y`). -""" -mid_point(M::AbstractManifold, p, q, ::Any) = mid_point(M, p, q) -mid_point!(M::AbstractManifold, y, p, q, ::Any) = mid_point!(M, y, p, q) - -function mid_point(M::Circle, p, q, x) - if distance(M, p, q) ≈ π - X = 0.5 * log(M, p, q) - Y = log(M, p, x) - return exp(M, p, (sign(X) == sign(Y) ? 1 : -1) * X) - end - return mid_point(M, p, q) -end - -function mid_point(M::Sphere, p, q, x) - if isapprox(M, p, -q) - X = log(M, p, x) / distance(M, p, x) * π - else - X = log(M, p, q) - end - return exp(M, p, 0.5 * X) -end -function mid_point!(M::Sphere, y, p, q, x) - if isapprox(M, p, -q) - X = log(M, p, x) / distance(M, p, x) * π - else - X = log(M, p, q) - end - y .= exp(M, p, 0.5 * X) - return y -end - -function prox_TV2(::Euclidean, λ, pointTuple::Tuple{T,T,T}, p::Int=1) where {T} - w = @SVector [1.0, -2.0, 1.0] - x = SVector(pointTuple) - if p == 1 # Example 3.2 in Bergmann, Laus, Steidl, Weinmann, 2014. - m = min.(Ref(λ), abs.(x .* w) / (dot(w, w))) - s = sign.(sum(x .* w)) - return x .- m .* s .* w - elseif p == 2 # Theorem 3.6 ibd. - t = λ * sum(x .* w) / (1 + λ * dot(w, w)) - return x .- t .* w - else - throw( - ErrorException( - "Proximal Map of TV2(Euclidean,λ,pT,p) not implemented for p=$(p) (requires p=1 or 2)", - ), - ) - end -end @doc raw""" reflect(M, f, x) reflect!(M, q, f, x) @@ -87,7 +19,8 @@ reflect!(M::AbstractManifold, q, pr::Function, x) = reflect!(M, q, pr(x), x) @doc raw""" reflect(M, p, x) reflect!(M, q, p, x) -reflect the point `x` from the manifold `M` at point `p`, i.e. + +Reflect the point `x` from the manifold `M` at point `p`, i.e. ````math \operatorname{refl}_p(x) = \exp_p(-\log_p x). diff --git a/src/functions/proximal_maps.jl b/src/functions/proximal_maps.jl index 1b9712aff1..897362b866 100644 --- a/src/functions/proximal_maps.jl +++ b/src/functions/proximal_maps.jl @@ -316,7 +316,7 @@ function prox_TV2( "Proximal Map of TV2(M, λ, x, p) not implemented for p=$(p) (requires p=1) on general manifolds.", ), ) - PowX = SVector(x) + PowX = [x...] PowM = PowerManifold(M, NestedPowerRepresentation(), 3) xR = PowX F(M, x) = 1 / 2 * distance(M, PowX, x)^2 + λ * costTV2(M, x) @@ -338,7 +338,7 @@ function prox_TV2!( "Proximal Map of TV2(M, λ, x, p) not implemented for p=$(p) (requires p=1) on general manifolds.", ), ) - PowX = SVector(x) + PowX = [x...] PowM = PowerManifold(M, NestedPowerRepresentation(), 3) copyto!(M, y, PowX) F(M, x) = 1 / 2 * distance(M, PowX, x)^2 + λ * costTV2(M, x) diff --git a/src/ext/LineSearchesExt.jl b/src/helpers/LineSearchesTypes.jl similarity index 60% rename from src/ext/LineSearchesExt.jl rename to src/helpers/LineSearchesTypes.jl index d8347c6433..59c3ecdc78 100644 --- a/src/ext/LineSearchesExt.jl +++ b/src/helpers/LineSearchesTypes.jl @@ -58,62 +58,7 @@ function LineSearchesStepsize( ) end -function (cs::LineSearchesStepsize)( - mp::AbstractManoptProblem, - s::AbstractManoptSolverState, - i::Int, - η=-get_gradient(s); - fp=get_cost(mp, get_iterate(s)), - kwargs..., -) - M = get_manifold(mp) - p = get_iterate(s) - X = get_gradient(s) - p_tmp = copy(M, p) - X_tmp = copy(M, p, X) - Y_tmp = copy(M, p, X) - f = get_cost_function(get_objective(mp)) - dphi_0 = real(inner(M, p, X, η)) - - # guess initial alpha - α0 = 1.0 - - # perform actual linesearch - - function ϕ(α) - retract!(M, p_tmp, p, η, α, cs.retraction_method) - return f(M, p_tmp) - end - function dϕ(α) - retract!(M, p_tmp, p, η, α, cs.retraction_method) - get_gradient!(mp, X_tmp, p_tmp) - vector_transport_to!(M, Y_tmp, p, η, p_tmp, cs.vector_transport_method) - return real(inner(M, p_tmp, Y_tmp, Y_tmp)) - end - function ϕdϕ(α) - # TODO: optimize? - retract!(M, p_tmp, p, η, α, cs.retraction_method) - get_gradient!(mp, X_tmp, p_tmp) - vector_transport_to!(M, Y_tmp, p, η, p_tmp, cs.vector_transport_method) - phi = f(M, p_tmp) - dphi = real(inner(M, p_tmp, Y_tmp, Y_tmp)) - return (phi, dphi) - end - - try - α, fp = cs.linesearch(ϕ, dϕ, ϕdϕ, α0, fp, dphi_0) - return α - catch ex - if isa(ex, LineSearches.LineSearchException) - # maybe indicate failure? - return zero(dphi_0) - else - rethrow(ex) - end - end -end - -function show(io::IO, cs::LineSearchesStepsize) +function Base.show(io::IO, cs::LineSearchesStepsize) return print( io, "LineSearchesStepsize($(cs.linesearch); retraction_method=$(cs.retraction_method), vector_transport_method=$(cs.vector_transport_method))", diff --git a/src/helpers/checks.jl b/src/helpers/checks.jl index 258e88d59f..b44738bb77 100644 --- a/src/helpers/checks.jl +++ b/src/helpers/checks.jl @@ -1,3 +1,18 @@ + +function plot_slope end + +""" + plot_slope(x, y; slope=2, line_base=0, a=0, b=2.0, i=1,j=length(x)) + +Plot the result from the error check functions, e.g. +[`check_gradient`](@ref), [`check_differential`](@ref), [`check_Hessian`](@ref) +on data `x,y` with two comparison lines + +1) `line_base` + t`slope` as the global slope the plot should have +2) `a` + `b*t` on the interval [`x[i]`, `x[j]`] for some (best fitting) comparison slope +""" +plot_slope(x, y) + """ prepare_check_result(log_range, errors, slope) diff --git a/src/plans/alternating_gradient_plan.jl b/src/plans/alternating_gradient_plan.jl index 4d1b4a8e2a..4d6c24c87c 100644 --- a/src/plans/alternating_gradient_plan.jl +++ b/src/plans/alternating_gradient_plan.jl @@ -11,7 +11,7 @@ An alternating gradient objective consists of !!! note - This Objective is usually defied using the `ProductManifold` from `Manifolds.jl`, so `Manifolds.jl` to be loaded. + This Objective is usually defined using the `ProductManifold` from `Manifolds.jl`, so `Manifolds.jl` to be loaded. # Constructors ManifoldAlternatingGradientObjective(F, gradF::Function; @@ -44,14 +44,6 @@ function ManifoldAlternatingGradientObjective( ) end -@doc raw""" - X = get_gradient(M::ProductManifold, ago::ManifoldAlternatingGradientObjective, p) - get_gradient!(M::ProductManifold, P::ManifoldAlternatingGradientObjective, X, p) - -Evaluate all summands gradients at a point `p` on the `ProductManifold M` (in place of `X`) -""" -get_gradient(M::ProductManifold, ::ManifoldAlternatingGradientObjective, ::Any...) - function get_gradient( M::AbstractManifold, mago::ManifoldAlternatingGradientObjective{AllocatingEvaluation,TC,<:Function}, @@ -59,13 +51,6 @@ function get_gradient( ) where {TC} return mago.gradient!!(M, p) end -function get_gradient( - M::AbstractManifold, - mago::ManifoldAlternatingGradientObjective{AllocatingEvaluation,TC,<:AbstractVector}, - p, -) where {TC} - return ProductRepr([gi(M, p) for gi in mago.gradient!!]...) -end function get_gradient!( M::AbstractManifold, X, @@ -111,32 +96,7 @@ function get_gradient!( mago.gradient!!(M, X, p) return X end -function get_gradient!( - M::AbstractManifold, - X, - mago::ManifoldAlternatingGradientObjective{InplaceEvaluation,TC,<:AbstractVector}, - p, -) where {TC} - for (gi, Xi) in zip(mago.gradient!!, submanifold_components(M, X)) - gi(M, Xi, p) - end - return X -end -@doc raw""" - X = get_gradient(M::AbstractManifold, p::ManifoldAlternatingGradientObjective, p, k) - get_gradient!(M::AbstractManifold, p::ManifoldAlternatingGradientObjective, X, p, k) - -Evaluate one of the component gradients ``\operatorname{grad}f_k``, ``k∈\{1,…,n\}``, at `x` (in place of `Y`). -""" -function get_gradient( - M::ProductManifold, - mago::ManifoldAlternatingGradientObjective{AllocatingEvaluation,TC,<:Function}, - p, - k, -) where {TC} - return get_gradient(M, mago, p)[M, k] -end function get_gradient( M::AbstractManifold, mago::ManifoldAlternatingGradientObjective{AllocatingEvaluation,TC,<:AbstractVector}, @@ -145,16 +105,6 @@ function get_gradient( ) where {TC} return mago.gradient!![k](M, p) end -function get_gradient!( - M::ProductManifold, - X, - mago::ManifoldAlternatingGradientObjective{AllocatingEvaluation,TC,<:Function}, - p, - k, -) where {TC} - copyto!(M[k], X, mago.gradient!!(M, p)[M, k]) - return X -end function get_gradient!( M::AbstractManifold, X, diff --git a/src/plans/cache.jl b/src/plans/cache.jl index 3372ed1f41..4d1b3efc79 100644 --- a/src/plans/cache.jl +++ b/src/plans/cache.jl @@ -1,8 +1,8 @@ # -# A Cache for Objectives +# A Simple Cache for Objectives # @doc raw""" - SimpleCacheObjective{O<:AbstractManifoldGradientObjective{E,TC,TG}, P, T,C} <: AbstractManifoldGradientObjective{E,TC,TG} + SimpleManifoldCachedObjective{O<:AbstractManifoldGradientObjective{E,TC,TG}, P, T,C} <: AbstractManifoldGradientObjective{E,TC,TG} Provide a simple cache for an [`AbstractManifoldGradientObjective`](@ref) that is for a given point `p` this cache stores a point `p` and a gradient ``\operatorname{grad} f(p)`` in `X` as well as a cost value ``f(p)`` in `c`. @@ -11,7 +11,7 @@ Both `X` and `c` are accompanied by booleans to keep track of their validity. # Constructor - SimpleCacheObjective(M::AbstractManifold, obj::AbstractManifoldGradientObjective; kwargs...) + SimpleManifoldCachedObjective(M::AbstractManifold, obj::AbstractManifoldGradientObjective; kwargs...) ## Keyword * `p` (`rand(M)`) – a point on the manifold to initialize the cache with @@ -19,9 +19,9 @@ Both `X` and `c` are accompanied by booleans to keep track of their validity. * `c` (`get_cost(M, obj, p)` or `0.0`) – a value to store the cost function in `initialize` * `initialized` (`true`) – whether to initialize the cached `X` and `c` or not. """ -mutable struct SimpleCacheObjective{ +mutable struct SimpleManifoldCachedObjective{ E<:AbstractEvaluationType,TC,TG,O<:AbstractManifoldGradientObjective{E,TC,TG},P,T,C -} <: AbstractManifoldGradientObjective{E,TC,TG} +} <: AbstractDecoratedManifoldObjective{E,O} objective::O p::P # a point X::T # a vector @@ -29,9 +29,8 @@ mutable struct SimpleCacheObjective{ c::C # a value for the cost c_valid::Bool end -dispatch_objective_decorator(::SimpleCacheObjective) = Val(true) -function SimpleCacheObjective( +function SimpleManifoldCachedObjective( M::AbstractManifold, obj::O; initialized=true, @@ -40,15 +39,12 @@ function SimpleCacheObjective( c=initialized ? get_cost(M, obj, p) : 0.0, ) where {E<:AbstractEvaluationType,TC,TG,O<:AbstractManifoldGradientObjective{E,TC,TG}} q = copy(M, p) - return SimpleCacheObjective{E,TC,TG,O,typeof(q),typeof(X),typeof(c)}( + return SimpleManifoldCachedObjective{E,TC,TG,O,typeof(q),typeof(X),typeof(c)}( obj, q, X, initialized, c, initialized ) end - -# -# Default implementation -# -function get_cost(M::AbstractManifold, sco::SimpleCacheObjective, p) +# Default implementations +function get_cost(M::AbstractManifold, sco::SimpleManifoldCachedObjective, p) scop_neq_p = sco.p != p if scop_neq_p || !sco.c_valid # else evaluate cost, invalidate grad if p changed @@ -60,9 +56,12 @@ function get_cost(M::AbstractManifold, sco::SimpleCacheObjective, p) end return sco.c end -get_cost_function(sco::SimpleCacheObjective) = get_cost_function(sco.objective) +get_cost_function(sco::SimpleManifoldCachedObjective) = (M, p) -> get_cost(M, sco, p) +function get_gradient_function(sco::SimpleManifoldCachedObjective) + return (M, p) -> get_gradient(M, sco, p) +end -function get_gradient(M::AbstractManifold, sco::SimpleCacheObjective, p) +function get_gradient(M::AbstractManifold, sco::SimpleManifoldCachedObjective, p) scop_neq_p = sco.p != p if scop_neq_p || !sco.X_valid get_gradient!(M, sco.X, sco.objective, p) @@ -73,7 +72,7 @@ function get_gradient(M::AbstractManifold, sco::SimpleCacheObjective, p) end return copy(M, p, sco.X) end -function get_gradient!(M::AbstractManifold, X, sco::SimpleCacheObjective, p) +function get_gradient!(M::AbstractManifold, X, sco::SimpleManifoldCachedObjective, p) scop_neq_p = sco.p != p if scop_neq_p || !sco.X_valid get_gradient!(M, sco.X, sco.objective, p) @@ -85,14 +84,15 @@ function get_gradient!(M::AbstractManifold, X, sco::SimpleCacheObjective, p) copyto!(M, X, sco.p, sco.X) return X end -get_gradient_function(sco::SimpleCacheObjective) = get_gradient_function(sco.objective) # # CostGradImplementation # function get_cost( M::AbstractManifold, - sco::SimpleCacheObjective{AllocatingEvaluation,TC,TG,<:ManifoldCostGradientObjective}, + sco::SimpleManifoldCachedObjective{ + AllocatingEvaluation,TC,TG,<:ManifoldCostGradientObjective + }, p, ) where {TC,TG} scop_neq_p = sco.p != p @@ -106,7 +106,9 @@ function get_cost( end function get_cost( M::AbstractManifold, - sco::SimpleCacheObjective{InplaceEvaluation,TC,TG,<:ManifoldCostGradientObjective}, + sco::SimpleManifoldCachedObjective{ + InplaceEvaluation,TC,TG,<:ManifoldCostGradientObjective + }, p, ) where {TC,TG} scop_neq_p = sco.p != p @@ -120,7 +122,9 @@ function get_cost( end function get_gradient( M::AbstractManifold, - sco::SimpleCacheObjective{AllocatingEvaluation,TC,TG,<:ManifoldCostGradientObjective}, + sco::SimpleManifoldCachedObjective{ + AllocatingEvaluation,TC,TG,<:ManifoldCostGradientObjective + }, p, ) where {TC,TG} scop_neq_p = sco.p != p @@ -135,7 +139,9 @@ function get_gradient( end function get_gradient( M::AbstractManifold, - sco::SimpleCacheObjective{InplaceEvaluation,TC,TG,<:ManifoldCostGradientObjective}, + sco::SimpleManifoldCachedObjective{ + InplaceEvaluation,TC,TG,<:ManifoldCostGradientObjective + }, p, ) where {TC,TG} scop_neq_p = sco.p != p @@ -151,7 +157,9 @@ end function get_gradient!( M::AbstractManifold, X, - sco::SimpleCacheObjective{AllocatingEvaluation,TC,TG,<:ManifoldCostGradientObjective}, + sco::SimpleManifoldCachedObjective{ + AllocatingEvaluation,TC,TG,<:ManifoldCostGradientObjective + }, p, ) where {TC,TG} scop_neq_p = sco.p != p @@ -167,7 +175,9 @@ end function get_gradient!( M::AbstractManifold, X, - sco::SimpleCacheObjective{InplaceEvaluation,TC,TG,<:ManifoldCostGradientObjective}, + sco::SimpleManifoldCachedObjective{ + InplaceEvaluation,TC,TG,<:ManifoldCostGradientObjective + }, p, ) where {TC,TG} scop_neq_p = sco.p != p @@ -181,6 +191,551 @@ function get_gradient!( return X end +# +# ManifoldCachedObjective constructor which errors by default since we can only define init +# for LRU Caches +# +@doc raw""" + ManifoldCachedObjective{E,P,O<:AbstractManifoldObjective{<:E},C<:NamedTuple{}} <: AbstractDecoratedManifoldObjective{E,P} + +Create a cache for an objective, based on a `NamedTuple` that stores some kind of cache. + +# Constructor + + ManifoldCachedObjective(M, o::AbstractManifoldObjective, caches::Vector{Symbol}; kwargs...) + +Create a cache for the [`AbstractManifoldObjective`](@ref) where the Symbols in `caches` indicate, +which function evaluations to cache. + +# Supported Symbols + +| Symbol | Caches calls to (incl. `!` variants) | Comment +| :-------------------------- | :------------------------------------- | :------------------------ | +| `:Constraints` | [`get_constraints`](@ref) | vector of numbers | +| `:Cost` | [`get_cost`](@ref) | | +| `:EqualityConstraint` | [`get_equality_constraint`](@ref) | numbers per (p,i) | +| `:EqualityConstraints` | [`get_equality_constraints`](@ref) | vector of numbers | +| `:GradEqualityConstraint` | [`get_grad_equality_constraint`](@ref) | tangent vector per (p,i) | +| `:GradEqualityConstraints` | [`get_grad_equality_constraints`](@ref)| vector of tangent vectors | +| `:GradInequalityConstraint` | [`get_inequality_constraint`](@ref) | tangent vector per (p,i) | +| `:GradInequalityConstraints`| [`get_inequality_constraints`](@ref) | vector of tangent vectors | +| `:Gradient` | [`get_gradient`](@ref)`(M,p)` | tangent vectors | +| `:Hessian` | [`get_hessian`](@ref) | tangent vectors | +| `:InequalityConstraint` | [`get_inequality_constraint`](@ref) | numbers per (p,j) | +| `:InequalityConstraints` | [`get_inequality_constraints`](@ref) | vector of numbers | +| `:Preconditioner` | [`get_preconditioner`](@ref) | tangent vectors | +| `:ProximalMap` | [`get_proximal_map`](@ref) | point per `(p,λ,i)` | +| `:StochasticGradients` | [`get_gradients`](@ref) | vector of tangent vectors | +| `:StochasticGradient` | [`get_gradient`](@ref)`(M, p, i)` | tangent vector per (p,i) | +| `:SubGradient` | [`get_subgradient`](@ref) | tangent vectors | +| `:SubtrahendGradient` | [`get_subtrahend_gradient`](@ref) | tangent vectors | + +# Keyword Arguments + +* `p` - (`rand(M)`) the type of the keys to be used in the caches. Defaults to the default representation on `M`. +* `value` - (`get_cost(M, objective, p)`) the type of values for numeric values in the cache, e.g. the cost +* `X` - (`zero_vector(M,p)`) the type of values to be cached for gradient and Hessian calls. +* `cache` - (`[:Cost]`) a vector of symbols indicating which function calls should be cached. +* `cache_size` - (`10`) number of (least recently used) calls to cache +* `cache_sizes` – (`Dict{Symbol,Int}()`) a named tuple or dictionary specifying the sizes individually for each cache. + + +""" +struct ManifoldCachedObjective{E,P,O<:AbstractManifoldObjective{<:E},C<:NamedTuple{}} <: + AbstractDecoratedManifoldObjective{E,P} + objective::O + cache::C +end +function ManifoldCachedObjective( + M::AbstractManifold, + objective::O, + caches::AbstractVector{<:Symbol}=[:Cost]; + p::P=rand(M), + value::R=get_cost(M, objective, p), + X::T=zero_vector(M, p), + cache_size::Int=10, + cache_sizes::Dict{Symbol,Int}=Dict{Symbol,Int}(), +) where {E,O<:AbstractManifoldObjective{E},R<:Real,P,T} + c = init_caches( + M, caches; p=p, value=value, X=X, cache_size=cache_size, cache_sizes=cache_sizes + ) + return ManifoldCachedObjective{E,O,O,typeof(c)}(objective, c) +end +function ManifoldCachedObjective( + M::AbstractManifold, + objective::O, + caches::AbstractVector{<:Symbol}=[:Cost]; + p::P=rand(M), + value::R=get_cost(M, objective, p), + X::T=zero_vector(M, p), + cache_size::Int=10, + cache_sizes::Dict{Symbol,Int}=Dict{Symbol,Int}(), +) where {E,O2,O<:AbstractDecoratedManifoldObjective{E,O2},R<:Real,P,T} + c = init_caches( + M, caches; p=p, value=value, X=X, cache_size=cache_size, cache_sizes=cache_sizes + ) + return ManifoldCachedObjective{E,O2,O,typeof(c)}(objective, c) +end + +""" + init_caches(M::AbstractManifold, caches, T; kwargs...) + +Given a vector of symbols `caches`, this function sets up the +`NamedTuple` of caches for points/vectors on `M`, +where `T` is the type of cache to use. +""" +function init_caches(M::AbstractManifold, caches, T=Nothing; kwargs...) + return throw(DomainError( + T, + """ + No function `init_caches` available for a caches $T. + For a good default load `LRUCache.jl`, for example: + Enter `using LRUCache`. + """, + )) +end + +# +# Default implementations - (a) check if field exists (b) try to get cache +function get_cost(M::AbstractManifold, co::ManifoldCachedObjective, p) + !(haskey(co.cache, :Cost)) && return get_cost(M, co.objective, p) + return get!(co.cache[:Cost], copy(M, p)) do + get_cost(M, co.objective, p) + end +end +get_cost_function(co::ManifoldCachedObjective) = (M, p) -> get_cost(M, co, p) +get_gradient_function(co::ManifoldCachedObjective) = (M, p) -> get_gradient(M, co, p) + +function get_gradient(M::AbstractManifold, co::ManifoldCachedObjective, p) + !(haskey(co.cache, :Gradient)) && return get_gradient(M, co.objective, p) + return copy( + M, + p, + get!(co.cache[:Gradient], copy(M, p)) do + get_gradient(M, co.objective, p) + end, + ) +end +function get_gradient!(M::AbstractManifold, X, co::ManifoldCachedObjective, p) + !(haskey(co.cache, :Gradient)) && return get_gradient!(M, X, co.objective, p) + copyto!( + M, + X, + p, + get!(co.cache[:Gradient], copy(M, p)) do + # This evaluates in place of X + get_gradient!(M, X, co.objective, p) + copy(M, p, X) #this creates a copy to be placed in the cache + end, #and we copy the values back to X + ) + return X +end + +# +# CostGradImplementation +function get_cost( + M::AbstractManifold, co::ManifoldCachedObjective{E,<:ManifoldCostGradientObjective}, p +) where {E<:AbstractEvaluationType} + #Neither cost not grad cached -> evaluate + all(.!(haskey.(Ref(co.cache), [:Cost, :Gradient]))) && + return get_cost(M, co.objective, p) + return get!(co.cache[:Cost], copy(M, p)) do + c, X = get_cost_and_gradient(M, co.objective, p) + #if this is evaluated, we can also set X + haskey(co.cache, :Gradient) && setindex!(co.cache[:Gradient], X, copy(M, p)) + c #but we also set the new cost here + end +end +function get_gradient( + M::AbstractManifold, co::ManifoldCachedObjective{E,<:ManifoldCostGradientObjective}, p +) where {E<:AllocatingEvaluation} + all(.!(haskey.(Ref(co.cache), [:Cost, :Gradient]))) && + return get_gradient(M, co.objective, p) + return copy( + M, + p, + get!(co.cache[:Gradient], copy(M, p)) do + c, X = get_cost_and_gradient(M, co.objective, p) + #if this is evaluated, we can also set c + haskey(co.cache, :Cost) && setindex!(co.cache[:Cost], c, copy(M, p)) + X #but we also set the new cost here + end, + ) +end +function get_gradient!( + M::AbstractManifold, + X, + co::ManifoldCachedObjective{E,<:ManifoldCostGradientObjective}, + p, +) where {E} + all(.!(haskey.(Ref(co.cache), [:Cost, :Gradient]))) && + return get_gradient!(M, X, co.objective, p) + return copyto!( + M, + X, + p, + get!(co.cache[:Gradient], copy(M, p)) do + c, _ = get_cost_and_gradient!(M, X, get_objective(co.objective), p) + #if this is evaluated, we can also set c + haskey(co.cache, :Cost) && setindex!(co.cache[:Cost], c, copy(M, p)) + X + end, + ) +end + +# +# Constraints +function get_constraints(M::AbstractManifold, co::ManifoldCachedObjective, p) + all(.!(haskey.(Ref(co.cache), [:Constraints]))) && + return get_constraints(M, co.objective, p) + return copy( + get!(co.cache[:Constraints], copy(M, p)) do + get_constraints(M, co.objective, p) + end, + ) +end +function get_equality_constraints(M::AbstractManifold, co::ManifoldCachedObjective, p) + all(.!(haskey.(Ref(co.cache), [:EqualityConstraints]))) && + return get_equality_constraints(M, co.objective, p) + return copy( + get!(co.cache[:EqualityConstraints], copy(M, p)) do + get_equality_constraints(M, co.objective, p) + end, + ) +end +function get_equality_constraint(M::AbstractManifold, co::ManifoldCachedObjective, p, i) + all(.!(haskey.(Ref(co.cache), [:EqualityConstraint]))) && + return get_equality_constraint(M, co.objective, p, i) + return copy( + get!(co.cache[:EqualityConstraint], (copy(M, p), i)) do + get_equality_constraint(M, co.objective, p, i) + end, + ) +end +function get_inequality_constraints(M::AbstractManifold, co::ManifoldCachedObjective, p) + all(.!(haskey.(Ref(co.cache), [:InequalityConstraints]))) && + return get_inequality_constraints(M, co.objective, p) + return copy( + get!(co.cache[:InequalityConstraints], copy(M, p)) do + get_inequality_constraints(M, co.objective, p) + end, + ) +end +function get_inequality_constraint(M::AbstractManifold, co::ManifoldCachedObjective, p, i) + all(.!(haskey.(Ref(co.cache), [:InequalityConstraint]))) && + return get_inequality_constraint(M, co.objective, p, i) + return copy( + get!(co.cache[:InequalityConstraint], (copy(M, p), i)) do + get_inequality_constraint(M, co.objective, p, i) + end, + ) +end +# +# Gradients of Constraints +function get_grad_equality_constraint( + M::AbstractManifold, co::ManifoldCachedObjective, p, j +) + !(haskey(co.cache, :GradEqualityConstraint)) && + return get_grad_equality_constraint(M, co.objective, p, j) + return copy( + M, + p, + get!(co.cache[:GradEqualityConstraint], (copy(M, p), j)) do + get_grad_equality_constraint(M, co.objective, p, j) + end, + ) +end +function get_grad_equality_constraint!( + M::AbstractManifold, X, co::ManifoldCachedObjective, p, j +) + !(haskey(co.cache, :GradEqualityConstraint)) && + return get_grad_equality_constraint!(M, X, co.objective, p, j) + copyto!( + M, + X, + p, + get!(co.cache[:GradEqualityConstraint], (copy(M, p), j)) do + # This evaluates in place of X + get_grad_equality_constraint!(M, X, co.objective, p, j) + copy(M, p, X) #this creates a copy to be placed in the cache + end, #and we copy the values back to X + ) + return X +end + +function get_grad_equality_constraints(M::AbstractManifold, co::ManifoldCachedObjective, p) + !(haskey(co.cache, :GradEqualityConstraints)) && + return get_grad_equality_constraints(M, co.objective, p) + return copy.( + Ref(M), + Ref(p), + get!(co.cache[:GradEqualityConstraints], copy(M, p)) do + get_grad_equality_constraints(M, co.objective, p) + end, + ) +end +function get_grad_equality_constraints!( + M::AbstractManifold, X, co::ManifoldCachedObjective, p +) + !(haskey(co.cache, :GradEqualityConstraints)) && + return get_grad_equality_constraints!(M, X, co.objective, p) + copyto!.( + Ref(M), + X, + Ref(p), + get!(co.cache[:GradEqualityConstraints], copy(M, p)) do + # This evaluates in place of X + get_grad_equality_constraints!(M, X, co.objective, p) + copy.(Ref(M), Ref(p), X) #this creates a copy to be placed in the cache + end, #and we copy the values back to X + ) + return X +end + +function get_grad_inequality_constraint( + M::AbstractManifold, co::ManifoldCachedObjective, p, j +) + !(haskey(co.cache, :GradInequalityConstraint)) && + return get_grad_inequality_constraint(M, co.objective, p, j) + return copy( + M, + p, + get!(co.cache[:GradInequalityConstraint], (copy(M, p), j)) do + get_grad_inequality_constraint(M, co.objective, p, j) + end, + ) +end +function get_grad_inequality_constraint!( + M::AbstractManifold, X, co::ManifoldCachedObjective, p, j +) + !(haskey(co.cache, :GradInequalityConstraint)) && + return get_grad_inequality_constraint!(M, X, co.objective, p, j) + copyto!( + M, + X, + p, + get!(co.cache[:GradInequalityConstraint], (copy(M, p), j)) do + # This evaluates in place of X + get_grad_inequality_constraint!(M, X, co.objective, p, j) + copy(M, p, X) #this creates a copy to be placed in the cache + end, #and we copy the values back to X + ) + return X +end + +function get_grad_inequality_constraints( + M::AbstractManifold, co::ManifoldCachedObjective, p +) + !(haskey(co.cache, :GradEqualityConstraints)) && + return get_grad_inequality_constraints(M, co.objective, p) + return copy.( + Ref(M), + Ref(p), + get!(co.cache[:GradInequalityConstraints], copy(M, p)) do + get_grad_inequality_constraints(M, co.objective, p) + end, + ) +end +function get_grad_inequality_constraints!( + M::AbstractManifold, X, co::ManifoldCachedObjective, p +) + !(haskey(co.cache, :GradInequalityConstraints)) && + return get_grad_inequality_constraints!(M, X, co.objective, p) + copyto!.( + Ref(M), + X, + Ref(p), + get!(co.cache[:GradInequalityConstraints], copy(M, p)) do + # This evaluates in place of X + get_grad_inequality_constraints!(M, X, co.objective, p) + copy.(Ref(M), Ref(p), X) #this creates a copy to be placed in the cache + end, #and we copy the values back to X + ) + return X +end + +# +# Hessian +function get_hessian(M::AbstractManifold, co::ManifoldCachedObjective, p, X) + !(haskey(co.cache, :Hessian)) && return get_hessian(M, co.objective, p, X) + return copy( + M, + p, + get!(co.cache[:Hessian], (copy(M, p), copy(M, p, X))) do + get_hessian(M, co.objective, p, X) + end, + ) +end +function get_hessian!(M::AbstractManifold, Y, co::ManifoldCachedObjective, p, X) + !(haskey(co.cache, :Hessian)) && return get_hessian!(M, Y, co.objective, p, X) + copyto!( + M, + Y, + p, #for the tricks performed here see get_gradient! + get!(co.cache[:Hessian], (copy(M, p), copy(M, p, X))) do + get_hessian!(M, Y, co.objective, p, X) + copy(M, p, Y) #store a copy of Y + end, + ) + return Y +end + +# +# Preconditioner +function get_preconditioner(M::AbstractManifold, co::ManifoldCachedObjective, p, X) + !(haskey(co.cache, :Preconditioner)) && return get_preconditioner(M, co.objective, p, X) + return copy( + M, + p, + get!(co.cache[:Preconditioner], (copy(M, p), copy(M, p, X))) do + get_preconditioner(M, co.objective, p, X) + end, + ) +end +function get_preconditioner!(M::AbstractManifold, Y, co::ManifoldCachedObjective, p, X) + !(haskey(co.cache, :Preconditioner)) && + return get_preconditioner!(M, Y, co.objective, p, X) + copyto!( + M, + Y, + p, #for the tricks performed here see get_gradient! + get!(co.cache[:Preconditioner], (copy(M, p), copy(M, p, X))) do + get_preconditioner!(M, Y, co.objective, p, X) + copy(M, p, Y) + end, + ) + return Y +end + +# +# Proximal Map +function get_proximal_map(M::AbstractManifold, co::ManifoldCachedObjective, λ, p, i) + !(haskey(co.cache, :ProximalMap)) && return get_proximal_map(M, co.objective, λ, p, i) + return copy( + M, + get!(co.cache[:ProximalMap], (copy(M, p), λ, i)) do #use the tuple (p,i) as key + get_proximal_map(M, co.objective, λ, p, i) + end, + ) +end +function get_proximal_map!(M::AbstractManifold, q, co::ManifoldCachedObjective, λ, p, i) + !(haskey(co.cache, :ProximalMap)) && + return get_proximal_map!(M, q, co.objective, λ, p, i) + copyto!( + M, + q, + get!(co.cache[:ProximalMap], (copy(M, p), λ, i)) do + get_proximal_map!(M, q, co.objective, λ, p, i) #compute inplace of q + copy(M, q) #store copy of q + end, + ) + return q +end +# +# Stochastic Gradient +function get_gradient(M::AbstractManifold, co::ManifoldCachedObjective, p, i) + !(haskey(co.cache, :StochasticGradient)) && return get_gradient(M, co.objective, p, i) + return copy( + M, + p, + get!(co.cache[:StochasticGradient], (copy(M, p), i)) do + get_gradient(M, co.objective, p, i) + end, + ) +end +function get_gradient!(M::AbstractManifold, X, co::ManifoldCachedObjective, p, i) + !(haskey(co.cache, :StochasticGradient)) && + return get_gradient!(M, X, co.objective, p, i) + copyto!( + M, + X, + p, + get!(co.cache[:StochasticGradient], (copy(M, p), i)) do + # This evaluates in place of X + get_gradient!(M, X, co.objective, p, i) + copy(M, p, X) #this creates a copy to be placed in the cache + end, #and we copy the values back to X + ) + return X +end + +function get_gradients(M::AbstractManifold, co::ManifoldCachedObjective, p) + !(haskey(co.cache, :StochasticGradients)) && return get_gradients(M, co.objective, p) + return copy.( + Ref(M), + Ref(p), + get!(co.cache[:StochasticGradients], copy(M, p)) do + get_gradients(M, co.objective, p) + end, + ) +end +function get_gradients!(M::AbstractManifold, X, co::ManifoldCachedObjective, p) + !(haskey(co.cache, :StochasticGradients)) && return get_gradients(M, X, co.objective, p) + copyto!.( + Ref(M), + X, + Ref(p), + get!(co.cache[:StochasticGradients], copy(M, p)) do + # This evaluates in place of X + get_gradients!(M, X, co.objective, p) + copy.(Ref(M), Ref(p), X) #this creates a copy to be placed in the cache + end, #and we copy the values back to X + ) + return X +end +# +# Subgradient +function get_subgradient(M::AbstractManifold, co::ManifoldCachedObjective, p) + !(haskey(co.cache, :SubGradient)) && return get_subgradient(M, co.objective, p) + return copy( + M, + p, + get!(co.cache[:SubGradient], copy(M, p)) do + get_subgradient(M, co.objective, p) + end, + ) +end +function get_subgradient!(M::AbstractManifold, X, co::ManifoldCachedObjective, p) + !(haskey(co.cache, :SubGradient)) && return get_subgradient!(M, X, co.objective, p) + copyto!( + M, + X, + p, #for the tricks performed here see get_gradient! + get!(co.cache[:SubGradient], copy(M, p)) do + get_subgradient!(M, X, co.objective, p) + copy(M, p, X) + end, + ) + return X +end + +# +# Subtrahend gradient (from DC) +function get_subtrahend_gradient(M::AbstractManifold, co::ManifoldCachedObjective, p) + !(haskey(co.cache, :SubtrahendGradient)) && + return get_subtrahend_gradient(M, co.objective, p) + return copy( + M, + p, + get!(co.cache[:SubtrahendGradient], copy(M, p)) do + get_subtrahend_gradient(M, co.objective, p) + end, + ) +end +function get_subtrahend_gradient!(M::AbstractManifold, X, co::ManifoldCachedObjective, p) + !(haskey(co.cache, :SubtrahendGradient)) && + return get_subtrahend_gradient!(M, X, co.objective, p) + copyto!( + M, + X, + p, #for the tricks performed here see get_gradient! + get!(co.cache[:SubtrahendGradient], copy(M, p)) do + get_subtrahend_gradient!(M, X, co.objective, p) + copy(M, p, X) + end, + ) + return X +end # # Factory # @@ -192,24 +747,47 @@ on the `AbstractManifold M` based on the symbol `cache`. The following caches are available -* `:Simple` generates a [`SimpleCacheObjective`](@ref) +* `:Simple` generates a [`SimpleManifoldCachedObjective`](@ref) +* `:LRU` generates a [`ManifoldCachedObjective`](@ref) where you should use the form + `(:LRU, [:Cost, :Gradient])` to specify what should be cached or + `(:LRU, [:Cost, :Gradient], 100)` to specify the cache size. + Here this variant defaults to `(:LRU, [:Cost, :Gradient], 100)`, + i.e. to cache up to 100 cost and gradient values.[^1] + +[^1]: + This cache requires [`LRUCache.jl`](https://github.com/JuliaCollections/LRUCache.jl) to be loaded as well. """ function objective_cache_factory(M, o, cache::Symbol) - (cache === :Simple) && return SimpleCacheObjective(M, o) + (cache === :Simple) && return SimpleManifoldCachedObjective(M, o) + (cache === :LRU) && + return ManifoldCachedObjective(M, o, [:Cost, :Gradient]; cache_size=100) return o end @doc raw""" - objective_cache_factory(M::AbstractManifold, o::AbstractManifoldObjective, cache::Tuple{Symbol, Array}Symbol) + objective_cache_factory(M::AbstractManifold, o::AbstractManifoldObjective, cache::Tuple{Symbol, Array, Array}) + objective_cache_factory(M::AbstractManifold, o::AbstractManifoldObjective, cache::Tuple{Symbol, Array}) Generate a cached variant of the [`AbstractManifoldObjective`](@ref) `o` on the `AbstractManifold M` based on the symbol `cache[1]`, -where the second element `cache[2]` is an array of further arguments for the cache and -the third is passed down as keyword arguments. +where the second element `cache[2]` are further arguments to the cache and +the optional third is passed down as keyword arguments. -For all availabel caches see the simpler variant with symbols. +For all available caches see the simpler variant with symbols. """ -function objective_cache_factory(M, o, cache::Tuple{Symbol,<:AbstractArray,<:AbstractArray}) - (cache[1] === :Simple) && return SimpleCacheObjective(M, o; cache[3]...) +function objective_cache_factory(M, o, cache::Tuple{Symbol,<:AbstractArray,I}) where {I} + (cache[1] === :Simple) && return SimpleManifoldCachedObjective(M, o; cache[3]...) + if (cache[1] === :LRU) + if (cache[3] isa Integer) + return ManifoldCachedObjective(M, o, cache[2]; cache_size=cache[3]) + else + return ManifoldCachedObjective(M, o, cache[2]; cache[3]...) + end + end + return o +end +function objective_cache_factory(M, o, cache::Tuple{Symbol,<:AbstractArray}) + (cache[1] === :Simple) && return SimpleManifoldCachedObjective(M, o) + (cache[1] === :LRU) && return ManifoldCachedObjective(M, o, cache[2]) return o end diff --git a/src/plans/constrained_plan.jl b/src/plans/constrained_plan.jl index f886bd2aa5..ca569a8285 100644 --- a/src/plans/constrained_plan.jl +++ b/src/plans/constrained_plan.jl @@ -271,6 +271,9 @@ containing the values of all constraints at `p`. function get_constraints(M::AbstractManifold, co::ConstrainedManifoldObjective, p) return [get_inequality_constraints(M, co, p), get_equality_constraints(M, co, p)] end +function get_constraints(M::AbstractManifold, admo::AbstractDecoratedManifoldObjective, p) + return get_constraints(M, get_objective(admo, false), p) +end function get_equality_constraints(mp::AbstractManoptProblem, p) return get_equality_constraints(get_manifold(mp), get_objective(mp), p) @@ -292,6 +295,11 @@ function get_equality_constraints( ) where {T<:AbstractEvaluationType} return [hj(M, p) for hj in co.h] end +function get_equality_constraints( + M::AbstractManifold, admo::AbstractDecoratedManifoldObjective, p +) + return get_equality_constraints(M, get_objective(admo, false), p) +end function get_equality_constraint(mp::AbstractManoptProblem, p, j) return get_equality_constraint(get_manifold(mp), get_objective(mp), p, j) @@ -315,6 +323,11 @@ function get_equality_constraint( ) where {T<:AbstractEvaluationType} return co.h[j](M, p) end +function get_equality_constraint( + M::AbstractManifold, admo::AbstractDecoratedManifoldObjective, p, j +) + return get_equality_constraint(M, get_objective(admo, false), p, j) +end function get_inequality_constraints(mp::AbstractManoptProblem, p) return get_inequality_constraints(get_manifold(mp), get_objective(mp), p) @@ -338,6 +351,12 @@ function get_inequality_constraints( return [gi(M, p) for gi in co.g] end +function get_inequality_constraints( + M::AbstractManifold, admo::AbstractDecoratedManifoldObjective, p +) + return get_inequality_constraints(M, get_objective(admo, false), p) +end + function get_inequality_constraint(mp::AbstractManoptProblem, p, i) return get_inequality_constraint(get_manifold(mp), get_objective(mp), p, i) end @@ -360,6 +379,11 @@ function get_inequality_constraint( ) where {T<:AbstractEvaluationType} return co.g[i](M, p) end +function get_inequality_constraint( + M::AbstractManifold, admo::AbstractDecoratedManifoldObjective, p, i +) + return get_inequality_constraint(M, get_objective(admo, false), p, i) +end function get_grad_equality_constraint(mp::AbstractManoptProblem, p, j) return get_grad_equality_constraint(get_manifold(mp), get_objective(mp), p, j) @@ -411,6 +435,11 @@ function get_grad_equality_constraint( co.grad_h!![j](M, X, p) return X end +function get_grad_equality_constraint( + M::AbstractManifold, admo::AbstractDecoratedManifoldObjective, p, j +) + return get_grad_equality_constraint(M, get_objective(admo, false), p, j) +end function get_grad_equality_constraint!(mp::AbstractManoptProblem, X, p, j) return get_grad_equality_constraint!(get_manifold(mp), X, get_objective(mp), p, j) @@ -470,6 +499,11 @@ function get_grad_equality_constraint!( co.grad_h!![j](M, X, p) return X end +function get_grad_equality_constraint!( + M::AbstractManifold, X, admo::AbstractDecoratedManifoldObjective, p, j +) + return get_grad_equality_constraint!(M, X, get_objective(admo, false), p, j) +end function get_grad_equality_constraints(mp::AbstractManoptProblem, p) return get_grad_equality_constraints(get_manifold(mp), get_objective(mp), p) @@ -518,6 +552,11 @@ function get_grad_equality_constraints( [grad_hi(M, Xj, p) for (Xj, grad_hi) in zip(X, co.grad_h!!)] return X end +function get_grad_equality_constraints( + M::AbstractManifold, admo::AbstractDecoratedManifoldObjective, p +) + return get_grad_equality_constraints(M, get_objective(admo, false), p) +end function get_grad_equality_constraints!(mp::AbstractManoptProblem, X, p) return get_grad_equality_constraints!(get_manifold(mp), X, get_objective(mp), p) @@ -568,6 +607,11 @@ function get_grad_equality_constraints!( end return X end +function get_grad_equality_constraints!( + M::AbstractManifold, X, admo::AbstractDecoratedManifoldObjective, p +) + return get_grad_equality_constraints!(M, X, get_objective(admo, false), p) +end function get_grad_inequality_constraint(mp::AbstractManoptProblem, p, i) return get_grad_inequality_constraint(get_manifold(mp), get_objective(mp), p, i) @@ -619,6 +663,11 @@ function get_grad_inequality_constraint( co.grad_g!![i](M, X, p) return X end +function get_grad_inequality_constraint( + M::AbstractManifold, admo::AbstractDecoratedManifoldObjective, p, i +) + return get_grad_inequality_constraint(M, get_objective(admo, false), p, i) +end function get_grad_inequality_constraint!(mp::AbstractManoptProblem, X, p, i) return get_grad_inequality_constraint!(get_manifold(mp), X, get_objective(mp), p, i) @@ -678,6 +727,11 @@ function get_grad_inequality_constraint!( co.grad_g!![i](M, X, p) return X end +function get_grad_inequality_constraint!( + M::AbstractManifold, X, admo::AbstractDecoratedManifoldObjective, p, i +) + return get_grad_inequality_constraint!(M, X, get_objective(admo, false), p, i) +end function get_grad_inequality_constraints(mp::AbstractManoptProblem, p) return get_grad_inequality_constraints(get_manifold(mp), get_objective(mp), p) @@ -726,6 +780,11 @@ function get_grad_inequality_constraints( [grad_gi(M, Xi, p) for (Xi, grad_gi) in zip(X, co.grad_g!!)] return X end +function get_grad_inequality_constraints( + M::AbstractManifold, admo::AbstractDecoratedManifoldObjective, p +) + return get_grad_inequality_constraints(M, get_objective(admo, false), p) +end function get_grad_inequality_constraints!(mp::AbstractManoptProblem, X, p) return get_grad_inequality_constraints!(get_manifold(mp), X, get_objective(mp), p) @@ -776,6 +835,11 @@ function get_grad_inequality_constraints!( end return X end +function get_grad_inequality_constraints!( + M::AbstractManifold, X, admo::AbstractDecoratedManifoldObjective, p +) + return get_grad_inequality_constraints!(M, X, get_objective(admo, false), p) +end function Base.show( io::IO, ::ConstrainedManifoldObjective{E,V} diff --git a/src/plans/cost_plan.jl b/src/plans/cost_plan.jl index f3c0e28fd9..f185dfb333 100644 --- a/src/plans/cost_plan.jl +++ b/src/plans/cost_plan.jl @@ -43,6 +43,9 @@ By default this implementation assumed that the cost is stored within `mco.cost` function get_cost(M::AbstractManifold, mco::AbstractManifoldCostObjective, p) return get_cost_function(mco)(M, p) end +function get_cost(M::AbstractManifold, admo::AbstractDecoratedManifoldObjective, p) + return get_cost(M, get_objective(admo, false), p) +end @doc raw""" get_cost_function(amco::AbstractManifoldCostObjective) @@ -50,3 +53,6 @@ end return the function to evaluate (just) the cost ``f(p)=c`` as a function `(M,p) -> c`. """ get_cost_function(mco::AbstractManifoldCostObjective) = mco.cost +function get_cost_function(admo::AbstractDecoratedManifoldObjective) + return get_cost_function(get_objective(admo, false)) +end diff --git a/src/plans/count.jl b/src/plans/count.jl new file mode 100644 index 0000000000..676c9c5923 --- /dev/null +++ b/src/plans/count.jl @@ -0,0 +1,428 @@ +""" + ManifoldCountObjective{E,P,O<:AbstractManifoldObjective,I<:Integer} <: AbstractDecoratedManifoldObjective{E,P} + +A wrapper for any [`AbstractManifoldObjective`](@ref) of type `O` to count different calls +to parts of the objective. + +# Fields + +* `counts` a dictionary of symbols mapping to integers keeping the counted values +* `objective` the wrapped objective + +# Supported Symbols + +| Symbol | Counts calls to (incl. `!` variants) | Comment | +| :-------------------------- | :------------------------------------- | :--------------------------- | +| `:Constraints` | [`get_constraints`](@ref) | | +| `:Cost` | [`get_cost`](@ref) | | +| `:EqualityConstraint` | [`get_equality_constraint`](@ref) | requires vector of counters | +| `:EqualityConstraints` | [`get_equality_constraints`](@ref) | does not count single access | +| `:GradEqualityConstraint` | [`get_grad_equality_constraint`](@ref) | requires vector of counters | +| `:GradEqualityConstraints` | [`get_grad_equality_constraints`](@ref)| does not count single access | +| `:GradInequalityConstraint` | [`get_inequality_constraint`](@ref) | requires vector of counters | +| `:GradInequalityConstraints`| [`get_inequality_constraints`](@ref) | does not count single access | +| `:Gradient` | [`get_gradient`](@ref)`(M,p)` | | +| `:Hessian` | [`get_hessian`](@ref) | | +| `:InequalityConstraint` | [`get_inequality_constraint`](@ref) | requires vector of counters | +| `:InequalityConstraints` | [`get_inequality_constraints`](@ref) | does not count single access | +| `:Preconditioner` | [`get_preconditioner`](@ref) | | +| `:ProximalMap` | [`get_proximal_map`](@ref) | | +| `:StochasticGradients` | [`get_gradients`](@ref) | | +| `:StochasticGradient` | [`get_gradient`](@ref)`(M, p, i)` | | +| `:SubGradient` | [`get_subgradient`](@ref) | | +| `:SubtrahendGradient` | [`get_subtrahend_gradient`](@ref) | | + +# Constructors + + ManifoldCountObjective(objective::AbstractManifoldObjective, counts::Dict{Symbol, <:Integer}) + +Initialise the `ManifoldCountObjective` to wrap `objective` initializing the set of counts + + ManifoldCountObjective(M::AbtractManifold, objective::AbstractManifoldObjective, count::AbstractVecor{Symbol}, init=0) + +Count function calls on `objective` using the symbols in `count` initialising all entries to `init`. +""" +struct ManifoldCountObjective{ + E,P,O<:AbstractManifoldObjective,I<:Union{<:Integer,AbstractVector{<:Integer}} +} <: AbstractDecoratedManifoldObjective{E,P} + counts::Dict{Symbol,I} + objective::O +end +function ManifoldCountObjective( + ::AbstractManifold, objective::O, counts::Dict{Symbol,I} +) where { + E<:AbstractEvaluationType, + I<:Union{<:Integer,AbstractVector{<:Integer}}, + O<:AbstractManifoldObjective{E}, +} + return ManifoldCountObjective{E,O,O,I}(counts, objective) +end +# Store the undecorated type of the input is decorated +function ManifoldCountObjective( + ::AbstractManifold, objective::O, counts::Dict{Symbol,I} +) where { + E<:AbstractEvaluationType, + I<:Union{<:Integer,AbstractVector{<:Integer}}, + P<:AbstractManifoldObjective, + O<:AbstractDecoratedManifoldObjective{E,P}, +} + return ManifoldCountObjective{E,P,O,I}(counts, objective) +end +function ManifoldCountObjective( + M::AbstractManifold, + objective::O, + count::AbstractVector{Symbol}, + init::I=0; + p::P=rand(M), +) where {P,I<:Integer,O<:AbstractManifoldObjective} + # Infere the sizes of the counters from the symbols if possible + counts = Pair{Symbol,Union{I,Vector{I}}}[] + for symbol in count + l = _get_counter_size(M, objective, symbol, p) + push!(counts, Pair(symbol, l == 1 ? init : fill(init, l))) + end + + return ManifoldCountObjective(M, objective, Dict(counts)) +end + +function _get_counter_size( + M::AbstractManifold, o::O, s::Symbol, p::P=rand(M) +) where {P,O<:AbstractManifoldObjective} + # vectorial counting cases + (s === :EqualityConstraint) && (return length(get_equality_constraints(M, o, p))) + (s === :GradEqualityConstraint) && (return length(get_equality_constraints(M, o, p))) + (s === :InequalityConstraint) && (return length(get_inequality_constraints(M, o, p))) + (s === :GradInequalityConstraint) && + (return length(get_inequality_constraints(M, o, p))) + # For now this only appears in ProximalMapObjective – so we can access its field + (s === :ProximalMap) && (return length(get_objective(o).proximal_maps!!)) + (s === :StochasticGradient) && (return length(get_gradients(M, o, p))) + return 1 #number - default +end + +function _count_if_exists(co::ManifoldCountObjective, s::Symbol) + return haskey(co.counts, s) && (co.counts[s] += 1) +end +function _count_if_exists(co::ManifoldCountObjective, s::Symbol, i) + if haskey(co.counts, s) + if (i == 1) && (ndims(co.counts[s]) == 0) + return co.counts[s] += 1 + elseif length(i) == ndims(co.counts[s]) && all(i .<= size(co.counts[s])) + return co.counts[s][i] += 1 + end + end +end + +""" + get_count(co::ManifoldCountObjective, s::Symbol, mode::Symbol=:None) + +Get the number of counts for a certain symbel `s`. + +Depending on the `mode` different results appear if the symbol does not exist in the dictionary + +* `:None` – (default) silent mode, returns `-1` for non-existing entries +* `:warn` – issues a warning if a field does not exist +* `:error` – issues an error if a field does not exist +""" +function get_count(co::ManifoldCountObjective, s::Symbol, mode::Symbol=:None) + if !haskey(co.counts, s) + msg = "There is no recorded count for $s." + (mode === :warn) && (@warn msg) + (mode === :error) && (error(msg)) + return -1 + end + return co.counts[s] +end +function get_count(o::AbstractManifoldObjective, s::Symbol, mode::Symbol=:None) + return _get_count(o, dispatch_objective_decorator(o), s, mode) +end +function _get_count(o::AbstractManifoldObjective, ::Val{false}, s, m) + return error("It seems $o does not provide access to a `ManifoldCountObjective`.") +end +function _get_count(o::AbstractManifoldObjective, ::Val{true}, s, m) + return get_count(get_objective(o, false), s, m) +end + +function get_count(co::ManifoldCountObjective, s::Symbol, i, mode::Symbol=:None) + if !haskey(co.counts, s) + msg = "There is no recorded count for :$s." + (mode === :warn) && (@warn msg) + (mode === :error) && (error(msg)) + return -1 + end + if !(ndims(i) == 0 && ndims(co.counts[s]) == 1) && ndims(i) != ndims(co.counts[s]) + msg = "The entry for :$s has $(ndims(co.counts[s])) dimensions but the index you provided has $(ndims(i))" + (mode === :warn) && (@warn msg) + (mode === :error) && (error(msg)) + return -1 + end + if ndims(i) == 0 && ndims(co.counts[s]) == 0 + if i > 1 + msg = "The entry for :$s is a number, but you provided the index $i > 1" + (mode === :warn) && (@warn msg) + (mode === :error) && (error(msg)) + return -1 + end + return co.counts[s] + end + if any(i .> size(co.counts[s])) + msg = "The index $i is out of range for the stored counts in :$s ($(size(co.counts[s])))." + (mode === :warn) && (@warn msg) + (mode === :error) && (error(msg)) + return -1 + end + return co.counts[s][i...] +end +function get_count(o::AbstractManifoldObjective, s::Symbol, i, mode::Symbol=:None) + return _get_count(o, dispatch_objective_decorator(o), s, i, mode) +end +function _get_count(o::AbstractManifoldObjective, ::Val{false}, s, i, m) + return error("It seems $o does not provide access to a `ManifoldCountObjective`.") +end +function _get_count(o::AbstractManifoldObjective, ::Val{true}, s, i, m) + return get_count(get_objective(o, false), s, i, m) +end + +""" + reset_counters(co::ManifoldCountObjective, value::Integer=0) + +Reset all values in the count objective to `value`. +""" +function reset_counters!(co::ManifoldCountObjective, value::Integer=0) + for s in keys(co.counts) + if (ndims(co.counts[s]) == 0) + co.counts[s] = value + else + co.counts[s] .= value + end + end + return co +end +function reset_counters!(o::AbstractDecoratedManifoldObjective, value::Integer=0) + return reset_counters!(get_objective(o, false), value) +end +function reset_counters!(o::AbstractManifoldObjective, value::Integer=0) + return error("It seems $o does not provide access to a `ManifoldCountObjective`.") +end + +# +# Overwrite accessors +# +function get_cost(M::AbstractManifold, co::ManifoldCountObjective, p) + _count_if_exists(co, :Cost) + return get_cost(M, co.objective, p) +end +function get_cost( + M::AbstractManifold, co::ManifoldCountObjective{E,<:ManifoldCostGradientObjective}, p +) where {E<:AbstractEvaluationType} + c, _ = get_cost_and_gradient(M, co, p) + return c +end +get_cost_function(co::ManifoldCountObjective) = (M, p) -> get_cost(M, co, p) + +function get_cost_and_gradient(M::AbstractManifold, co::ManifoldCountObjective, p) + _count_if_exists(co, :Cost) + _count_if_exists(co, :Gradient) + return get_cost_and_gradient(M, co.objective, p) +end + +function get_cost_and_gradient!(M::AbstractManifold, X, co::ManifoldCountObjective, p) + _count_if_exists(co, :Cost) + _count_if_exists(co, :Gradient) + return get_cost_and_gradient!(M, X, co.objective, p) +end + +get_gradient_function(co::ManifoldCountObjective) = (M, p) -> get_gradient(M, co, p) +function get_gradient(M::AbstractManifold, co::ManifoldCountObjective, p) + _count_if_exists(co, :Gradient) + return get_gradient(M, co.objective, p) +end +function get_gradient!(M::AbstractManifold, X, co::ManifoldCountObjective, p) + _count_if_exists(co, :Gradient) + get_gradient!(M, X, co.objective, p) + return X +end +function get_gradient( + M::AbstractManifold, co::ManifoldCountObjective{E,<:ManifoldCostGradientObjective}, p +) where {E<:AbstractEvaluationType} + _, X = get_cost_and_gradient(M, co, p) + return X +end +function get_gradient!( + M::AbstractManifold, X, co::ManifoldCountObjective{E,<:ManifoldCostGradientObjective}, p +) where {E<:AbstractEvaluationType} + get_cost_and_gradient!(M, X, co, p) + return X +end + +function get_hessian(M::AbstractManifold, co::ManifoldCountObjective, p, X) + _count_if_exists(co, :Hessian) + return get_hessian(M, co.objective, p, X) +end +function get_hessian!(M::AbstractManifold, Y, co::ManifoldCountObjective, p, X) + _count_if_exists(co, :Hessian) + get_hessian!(M, Y, co.objective, p, X) + return Y +end + +function get_preconditioner(M::AbstractManifold, co::ManifoldCountObjective, p, X) + _count_if_exists(co, :Preconditioner) + return get_preconditioner(M, co.objective, p, X) +end +function get_preconditioner!(M::AbstractManifold, Y, co::ManifoldCountObjective, p, X) + _count_if_exists(co, :Preconditioner) + get_preconditioner!(M, Y, co.objective, p, X) + return Y +end + +# +# Constraint +function get_constraints(M::AbstractManifold, co::ManifoldCountObjective, p) + _count_if_exists(co, :Constraints) + return get_constraints(M, co.objective, p) +end +function get_equality_constraints(M::AbstractManifold, co::ManifoldCountObjective, p) + _count_if_exists(co, :EqualityConstraints) + return get_equality_constraints(M, co.objective, p) +end +function get_equality_constraint(M::AbstractManifold, co::ManifoldCountObjective, p, i) + _count_if_exists(co, :EqualityConstraint, i) + return get_equality_constraint(M, co.objective, p, i) +end +function get_inequality_constraints(M::AbstractManifold, co::ManifoldCountObjective, p) + _count_if_exists(co, :InequalityConstraints) + return get_inequality_constraints(M, co.objective, p) +end +function get_inequality_constraint(M::AbstractManifold, co::ManifoldCountObjective, p, i) + _count_if_exists(co, :InequalityConstraint, i) + return get_inequality_constraint(M, co.objective, p, i) +end + +function get_grad_equality_constraints(M::AbstractManifold, co::ManifoldCountObjective, p) + _count_if_exists(co, :GradEqualityConstraints) + return get_grad_equality_constraints(M, co.objective, p) +end +function get_grad_equality_constraints!( + M::AbstractManifold, X, co::ManifoldCountObjective, p +) + _count_if_exists(co, :GradEqualityConstraints) + return get_grad_equality_constraints!(M, X, co.objective, p) +end +function get_grad_equality_constraint(M::AbstractManifold, co::ManifoldCountObjective, p, i) + _count_if_exists(co, :GradEqualityConstraint, i) + return get_grad_equality_constraint(M, co.objective, p, i) +end +function get_grad_equality_constraint!( + M::AbstractManifold, X, co::ManifoldCountObjective, p, i +) + _count_if_exists(co, :GradEqualityConstraint, i) + return get_grad_equality_constraint!(M, X, co.objective, p, i) +end +function get_grad_inequality_constraints(M::AbstractManifold, co::ManifoldCountObjective, p) + _count_if_exists(co, :GradInequalityConstraints) + return get_grad_inequality_constraints(M, co.objective, p) +end +function get_grad_inequality_constraints!( + M::AbstractManifold, X, co::ManifoldCountObjective, p +) + _count_if_exists(co, :GradInequalityConstraints) + return get_grad_inequality_constraints!(M, X, co.objective, p) +end +function get_grad_inequality_constraint( + M::AbstractManifold, co::ManifoldCountObjective, p, i +) + _count_if_exists(co, :GradInequalityConstraint, i) + return get_grad_inequality_constraint(M, co.objective, p, i) +end +function get_grad_inequality_constraint!( + M::AbstractManifold, X, co::ManifoldCountObjective, p, i +) + _count_if_exists(co, :GradInequalityConstraint, i) + return get_grad_inequality_constraint!(M, X, co.objective, p, i) +end + +# +# proxes +function get_proximal_map(M::AbstractManifold, co::ManifoldCountObjective, λ, p) + _count_if_exists(co, :ProximalMap) + return get_proximal_map(M, co.objective, λ, p) +end +function get_proximal_map!(M::AbstractManifold, q, co::ManifoldCountObjective, λ, p) + _count_if_exists(co, :ProximalMap) + return get_proximal_map!(M, q, co.objective, λ, p) +end +function get_proximal_map(M::AbstractManifold, co::ManifoldCountObjective, λ, p, i) + _count_if_exists(co, :ProximalMap, i) + return get_proximal_map(M, co.objective, λ, p, i) +end +function get_proximal_map!(M::AbstractManifold, q, co::ManifoldCountObjective, λ, p, i) + _count_if_exists(co, :ProximalMap, i) + return get_proximal_map!(M, q, co.objective, λ, p, i) +end + +# +# DC +function get_subtrahend_gradient(M::AbstractManifold, co::ManifoldCountObjective, p) + _count_if_exists(co, :SubtrahendGradient) + return get_subtrahend_gradient(M, co.objective, p) +end +function get_subtrahend_gradient!(M::AbstractManifold, X, co::ManifoldCountObjective, p) + _count_if_exists(co, :SubtrahendGradient) + return get_subtrahend_gradient!(M, X, co.objective, p) +end + +# +# Subgradient +function get_subgradient(M::AbstractManifold, co::ManifoldCountObjective, p) + _count_if_exists(co, :SubGradient) + return get_subgradient(M, co.objective, p) +end +function get_subgradient!(M::AbstractManifold, X, co::ManifoldCountObjective, p) + _count_if_exists(co, :SubGradient) + return get_subgradient!(M, X, co.objective, p) +end + +# +# Stochastic Gradient +function get_gradients(M::AbstractManifold, co::ManifoldCountObjective, p) + _count_if_exists(co, :StochasticGradients) + return get_gradients(M, co.objective, p) +end +function get_gradients!(M::AbstractManifold, X, co::ManifoldCountObjective, p) + _count_if_exists(co, :StochasticGradients) + return get_gradients!(M, X, co.objective, p) +end +function get_gradient(M::AbstractManifold, co::ManifoldCountObjective, p, i) + _count_if_exists(co, :StochasticGradient, i) + return get_gradient(M, co.objective, p, i) +end +function get_gradient!(M::AbstractManifold, X, co::ManifoldCountObjective, p, i) + _count_if_exists(co, :StochasticGradient, i) + return get_gradient!(M, X, co.objective, p, i) +end + +function objective_count_factory( + M::AbstractManifold, o::AbstractManifoldCostObjective, counts::Vector{<:Symbol} +) + return ManifoldCountObjective(M, o, counts) +end + +function status_summary(co::ManifoldCountObjective) + longest_key_length = max(length.(["$c" for c in keys(co.counts)])...) + s = "## Statistics on function calls\n" + count_strings = [ + " * :$(rpad("$(c[1])",longest_key_length)) : $(c[2])" for c in co.counts + ] + s2 = status_summary(co.objective) + !(co.objective isa AbstractDecoratedManifoldObjective) && (s2 = "on a $(s2)") + return "$(s)$(join(count_strings,"\n"))\n$s2" +end + +function show(io::IO, co::ManifoldCountObjective) + return print(io, "$(status_summary(co))") +end +function show( + io::IO, t::Tuple{<:ManifoldCountObjective,S} +) where {S<:AbstractManoptSolverState} + return print(io, "$(t[2])\n\n$(t[1])") +end diff --git a/src/plans/difference_of_convex_plan.jl b/src/plans/difference_of_convex_plan.jl index 6424c6e12d..ba842487fe 100644 --- a/src/plans/difference_of_convex_plan.jl +++ b/src/plans/difference_of_convex_plan.jl @@ -67,6 +67,12 @@ function get_subtrahend_gradient( X = zero_vector(M, p) return doco.∂h!!(M, X, p) end +function get_subtrahend_gradient( + M::AbstractManifold, admo::AbstractDecoratedManifoldObjective, p +) + return get_subtrahend_gradient(M, get_objective(admo, false), p) +end + function get_subtrahend_gradient!( M::AbstractManifold, X, @@ -80,6 +86,11 @@ function get_subtrahend_gradient!( ) return doco.∂h!!(M, X, p) end +function get_subtrahend_gradient!( + M::AbstractManifold, X, admo::AbstractDecoratedManifoldObjective, p +) + return get_subtrahend_gradient!(M, X, get_objective(admo, false), p) +end @doc raw""" LinearizedDCCost diff --git a/src/plans/gradient_plan.jl b/src/plans/gradient_plan.jl index d5a61d2013..bd250f533c 100644 --- a/src/plans/gradient_plan.jl +++ b/src/plans/gradient_plan.jl @@ -17,6 +17,9 @@ Depending on the [`AbstractEvaluationType`](@ref) `E` this is a function * `(M, X, p) -> X` for the [`InplaceEvaluation`](@ref), i.e. working inplace of `X`. """ get_gradient_function(amgo::AbstractManifoldGradientObjective) = amgo.gradient!! +function get_gradient_function(admo::AbstractDecoratedManifoldObjective) + return get_gradient_function(get_objective(admo, false)) +end @doc raw""" ManifoldGradientObjective{T<:AbstractEvaluationType} <: AbstractManifoldGradientObjective{T} @@ -87,8 +90,36 @@ function get_gradient_function(cgo::ManifoldCostGradientObjective) return (M, p) -> get_gradient(M, cgo, p) end +# +# and indernal helper to make the dispatch nicer +# +function get_cost_and_gradient( + M::AbstractManifold, cgo::ManifoldCostGradientObjective{AllocatingEvaluation}, p +) + return cgo.costgrad!!(M, p) +end +function get_cost_and_gradient( + M::AbstractManifold, cgo::ManifoldCostGradientObjective{InplaceEvaluation}, p +) + X = zero_vector(M, p) + return cgo.costgrad!!(M, X, p) +end + +function get_cost_and_gradient!( + M::AbstractManifold, X, cgo::ManifoldCostGradientObjective{AllocatingEvaluation}, p +) + (c, Y) = cgo.costgrad!!(M, p) + copyto!(M, X, p, Y) + return (c, X) +end +function get_cost_and_gradient!( + M::AbstractManifold, X, cgo::ManifoldCostGradientObjective{InplaceEvaluation}, p +) + return cgo.costgrad!!(M, X, p) +end + function get_cost(M::AbstractManifold, cgo::ManifoldCostGradientObjective, p) - v, _ = cgo.costgrad!!(M, p) + v, _ = get_cost_and_gradient(M, cgo, p) return v end @@ -124,6 +155,9 @@ vector `X` comes second. """ get_gradient(M::AbstractManifold, mgo::AbstractManifoldGradientObjective, p) +function get_gradient(M::AbstractManifold, admo::AbstractDecoratedManifoldObjective, p) + return get_gradient(M, get_objective(admo, false), p) +end function get_gradient( M::AbstractManifold, mgo::AbstractManifoldGradientObjective{AllocatingEvaluation}, p ) @@ -136,18 +170,13 @@ function get_gradient( mgo.gradient!!(M, X, p) return X end -function get_gradient( - M::AbstractManifold, mgo::ManifoldCostGradientObjective{AllocatingEvaluation}, p -) - _, X = mgo.costgrad!!(M, p) +function get_gradient(M::AbstractManifold, mcgo::ManifoldCostGradientObjective, p) + _, X = get_cost_and_gradient(M, mcgo, p) return X end -function get_gradient( - M::AbstractManifold, mgo::ManifoldCostGradientObjective{InplaceEvaluation}, p -) - X = zero_vector(M, p) - mgo.costgrad!!(M, X, p) - return X + +function get_gradient!(M::AbstractManifold, X, admo::AbstractDecoratedManifoldObjective, p) + return get_gradient!(M, X, get_objective(admo, false), p) end function get_gradient!( @@ -162,17 +191,8 @@ function get_gradient!( mgo.gradient!!(M, X, p) return X end -function get_gradient!( - M::AbstractManifold, X, mgo::ManifoldCostGradientObjective{AllocatingEvaluation}, p -) - _, Y = mgo.costgrad!!(M, p) - copyto!(M, X, p, Y) - return X -end -function get_gradient!( - M::AbstractManifold, X, mgo::ManifoldCostGradientObjective{InplaceEvaluation}, p -) - mgo.costgrad!!(M, X, p) +function get_gradient!(M::AbstractManifold, X, mcgo::ManifoldCostGradientObjective, p) + get_cost_and_gradient!(M, X, mcgo, p) return X end diff --git a/src/plans/hessian_plan.jl b/src/plans/hessian_plan.jl index 3a300c80f9..2fed2ebde5 100644 --- a/src/plans/hessian_plan.jl +++ b/src/plans/hessian_plan.jl @@ -73,6 +73,10 @@ end function get_hessian!(amp::AbstractManoptProblem, Y, p, X) return get_hessian!(get_manifold(amp), Y, get_objective(amp), p, X) end + +function get_hessian(M::AbstractManifold, admo::AbstractDecoratedManifoldObjective, p, X) + return get_hessian(M, get_objective(admo, false), p, X) +end function get_hessian( M::AbstractManifold, mho::ManifoldHessianObjective{AllocatingEvaluation}, p, X ) @@ -85,6 +89,11 @@ function get_hessian( mho.hessian!!(M, Y, p, X) return Y end +function get_hessian!( + M::AbstractManifold, Y, admo::AbstractDecoratedManifoldObjective, p, X +) + return get_hessian!(M, Y, get_objective(admo, false), p, X) +end function get_hessian!( M::AbstractManifold, Y, mho::ManifoldHessianObjective{AllocatingEvaluation}, p, X ) @@ -133,12 +142,23 @@ function get_preconditioner( mho.preconditioner!!(M, Y, p, X) return Y end +function get_preconditioner( + M::AbstractManifold, admo::AbstractDecoratedManifoldObjective, p, X +) + return get_preconditioner(M, get_objective(admo, false), p, X) +end + function get_preconditioner!( M::AbstractManifold, Y, mho::ManifoldHessianObjective{AllocatingEvaluation}, p, X ) copyto!(M, Y, p, mho.preconditioner!!(M, p, X)) return Y end +function get_preconditioner!( + M::AbstractManifold, Y, admo::AbstractDecoratedManifoldObjective, p, X +) + return get_preconditioner!(M, Y, get_objective(admo, false), p, X) +end function get_preconditioner!( M::AbstractManifold, Y, mho::ManifoldHessianObjective{InplaceEvaluation}, p, X ) diff --git a/src/plans/higher_order_primal_dual_plan.jl b/src/plans/higher_order_primal_dual_plan.jl index f1f32aba55..9dc17596fe 100644 --- a/src/plans/higher_order_primal_dual_plan.jl +++ b/src/plans/higher_order_primal_dual_plan.jl @@ -238,6 +238,11 @@ function get_differential_primal_prox( pdsno.diff_prox_f!!(M, Y, σ, p, X) return Y end +function get_differential_primal_prox( + M::AbstractManifold, admo::AbstractDecoratedManifoldObjective, σ, p, X +) + return get_differential_primal_prox(M, get_objective(admo, false), σ, p, X) +end function get_differential_primal_prox!( M::AbstractManifold, Y, @@ -260,6 +265,11 @@ function get_differential_primal_prox!( pdsno.diff_prox_f!!(M, Y, σ, p, X) return Y end +function get_differential_primal_prox!( + M::AbstractManifold, Y, admo::AbstractDecoratedManifoldObjective, σ, p, X +) + return get_differential_primal_prox!(M, Y, get_objective(admo, false), σ, p, X) +end @doc raw""" η = get_differential_dual_prox(N::AbstractManifold, pdsno::PrimalDualManifoldSemismoothNewtonObjective, n, τ, X, ξ) @@ -313,6 +323,11 @@ function get_differential_dual_prox( pdsno.diff_prox_g_dual!!(N, η, n, τ, X, ξ) return η end +function get_differential_dual_prox( + M::AbstractManifold, admo::AbstractDecoratedManifoldObjective, n, τ, X, ξ +) + return get_differential_dual_prox(M, get_objective(admo, false), n, τ, X, ξ) +end function get_differential_dual_prox!( N::AbstractManifold, η, @@ -337,3 +352,8 @@ function get_differential_dual_prox!( pdsno.diff_prox_g_dual!!(N, η, n, τ, X, ξ) return η end +function get_differential_dual_prox!( + M::AbstractManifold, η, admo::AbstractDecoratedManifoldObjective, n, τ, X, ξ +) + return get_differential_dual_prox!(M, η, get_objective(admo, false), n, τ, X, ξ) +end diff --git a/src/plans/objective.jl b/src/plans/objective.jl index 32d8cac722..d56a32aa3b 100644 --- a/src/plans/objective.jl +++ b/src/plans/objective.jl @@ -6,7 +6,7 @@ An abstract type to specify the kind of evaluation a [`AbstractManifoldObjective abstract type AbstractEvaluationType end @doc raw""" - AbstractManifoldObjective{T<:AbstractEvaluationType} + AbstractManifoldObjective{E<:AbstractEvaluationType} Describe the collection of the optimization function ``f\colon \mathcal M → \bbR` (or even a vectorial range) and its corresponding elements, which might for example be a gradient or (one or more) prxomial maps. @@ -20,7 +20,16 @@ All these elements should usually be implemented as functions the type `T` indicates the global [`AbstractEvaluationType`](@ref). """ -abstract type AbstractManifoldObjective{T<:AbstractEvaluationType} end +abstract type AbstractManifoldObjective{E<:AbstractEvaluationType} end + +@doc raw""" + AbstractDecoratedManifoldObjective{E<:AbstractEvaluationType,O<:AbstractManifoldObjective} + +A common supertype for all decorators of [`AbstractManifoldObjective`](@ref)s to simplify dispatch. + The second parameter should refer to the undecorated objective (i.e. the most inner one). +""" +abstract type AbstractDecoratedManifoldObjective{E,O<:AbstractManifoldObjective} <: + AbstractManifoldObjective{E} end @doc raw""" AllocatingEvaluation <: AbstractEvaluationType @@ -38,6 +47,25 @@ do not allocate memory but work on their input, i.e. in place. """ struct InplaceEvaluation <: AbstractEvaluationType end +struct ReturnManifoldObjective{E,P,O<:AbstractManifoldObjective{E}} <: + AbstractDecoratedManifoldObjective{E,P} + objective::O +end +function ReturnManifoldObjective( + o::O +) where {E<:AbstractEvaluationType,O<:AbstractManifoldObjective{E}} + return ReturnManifoldObjective{E,O,O}(o) +end +function ReturnManifoldObjective( + o::O +) where { + E<:AbstractEvaluationType, + P<:AbstractManifoldObjective, + O<:AbstractDecoratedManifoldObjective{E,P}, +} + return ReturnManifoldObjective{E,P,O}(o) +end + """ dispatch_objective_decorator(o::AbstractManoptSolverState) @@ -49,6 +77,7 @@ Decorators indicate this by returning `Val{true}` for further dispatch. The default is `Val{false}`, i.e. by default an state is not decorated. """ dispatch_objective_decorator(::AbstractManifoldObjective) = Val(false) +dispatch_objective_decorator(::AbstractDecoratedManifoldObjective) = Val(true) """ is_object_decorator(s::AbstractManifoldObjective) @@ -60,15 +89,48 @@ function is_objective_decorator(s::AbstractManifoldObjective) end @doc raw""" - get_objective(o::AbstractManifoldObjective) + get_objective(o::AbstractManifoldObjective, recursive=true) -return the undecorated [`AbstractManifoldObjective`](@ref) of the (possibly) decorated `o`. +return the (one step) undecorated [`AbstractManifoldObjective`](@ref) of the (possibly) decorated `o`. As long as your decorated objective stores the objetive within `o.objective` and the [`dispatch_objective_decorator`](@ref) is set to `Val{true}`, the internal state are extracted automatically. + +By default the objective that is stored within a decorated objective is assumed to be at +`o.objective`. Overwrtie `_get_objective(o, ::Val{true}, recursive) to change this bevahiour for your objective `o` +for both the recursive and the nonrecursive case. + +If `recursive` is set to `false`, only the most outer decorator is taken away instead of all. """ -function get_objective(o::AbstractManifoldObjective) - return _get_objective(o, dispatch_objective_decorator(o)) +function get_objective(o::AbstractManifoldObjective, recursive=true) + return _get_objective(o, dispatch_objective_decorator(o), recursive) +end +_get_objective(o::AbstractManifoldObjective, ::Val{false}, rec) = o +function _get_objective(o::AbstractManifoldObjective, ::Val{true}, rec) + return rec ? get_objective(o.objective) : o.objective +end +function status_summary(o::AbstractManifoldObjective{E}) where {E} + return "$(nameof(typeof(o))){$E}" +end +# Default undecorate for summary +function status_summary(co::AbstractDecoratedManifoldObjective) + return status_summary(get_objective(co, false)) +end + +function show(io::IO, o::AbstractManifoldObjective{E}) where {E} + return print(io, "$(nameof(typeof(o))){$E}") +end +# Default undecorate for show +function show(io::IO, co::AbstractDecoratedManifoldObjective) + return show(io, get_objective(co, false)) +end + +function show(io::IO, t::Tuple{<:AbstractManifoldObjective,P}) where {P} + return print( + io, + """ +$(status_summary(t[1])) + +To access the solver result, call `get_solver_result` on this variable.""", + ) end -_get_objective(o::AbstractManifoldObjective, ::Val{false}) = o -_get_objective(o::AbstractManifoldObjective, ::Val{true}) = get_objective(o.objective) diff --git a/src/plans/plan.jl b/src/plans/plan.jl index 5a2619b391..7700abde96 100644 --- a/src/plans/plan.jl +++ b/src/plans/plan.jl @@ -21,6 +21,7 @@ include("stepsize.jl") include("cost_plan.jl") include("gradient_plan.jl") +include("alternating_gradient_plan.jl") include("constrained_plan.jl") include("augmented_lagrangian_plan.jl") include("conjugate_gradient_plan.jl") @@ -41,3 +42,4 @@ include("stochastic_gradient_plan.jl") include("subsolver_plan.jl") include("cache.jl") +include("count.jl") diff --git a/src/plans/primal_dual_plan.jl b/src/plans/primal_dual_plan.jl index 691b852d8e..4e50f08c5a 100644 --- a/src/plans/primal_dual_plan.jl +++ b/src/plans/primal_dual_plan.jl @@ -144,6 +144,12 @@ function get_primal_prox( q = allocate_result(M, get_primal_prox, p) return apdmo.prox_f!!(M, q, σ, p) end +function get_primal_prox( + M::AbstractManifold, admo::AbstractDecoratedManifoldObjective, σ, p +) + return get_primal_prox(M, get_objective(admo, false), σ, p) +end + function get_primal_prox!( M::AbstractManifold, q, @@ -164,6 +170,11 @@ function get_primal_prox!( apdmo.prox_f!!(M, q, σ, p) return q end +function get_primal_prox!( + M::AbstractManifold, q, admo::AbstractDecoratedManifoldObjective, σ, p +) + return get_primal_prox!(M, q, get_objective(admo, false), σ, p) +end @doc raw""" Y = get_dual_prox(N::AbstractManifold, apdmo::AbstractPrimalDualManifoldObjective, n, τ, X) @@ -207,6 +218,12 @@ function get_dual_prox( apdmo.prox_g_dual!!(M, Y, n, τ, X) return Y end +function get_dual_prox( + M::AbstractManifold, admo::AbstractDecoratedManifoldObjective, n, τ, X +) + return get_dual_prox(M, get_objective(admo, false), n, τ, X) +end + function get_dual_prox!( M::AbstractManifold, Y, @@ -229,6 +246,11 @@ function get_dual_prox!( apdmo.prox_g_dual!!(M, Y, n, τ, X) return Y end +function get_dual_prox!( + M::AbstractManifold, Y, admo::AbstractDecoratedManifoldObjective, n, τ, X +) + return get_dual_prox!(M, Y, get_objective(admo, false), n, τ, X) +end @doc raw""" Y = linearized_forward_operator(M::AbstractManifold, N::AbstractManifold, apdmo::AbstractPrimalDualManifoldObjective, m, X, n) @@ -275,6 +297,17 @@ function linearized_forward_operator( apdmo.linearized_forward_operator!!(M, Y, m, X) return Y end +function linearized_forward_operator( + M::AbstractManifold, + N::AbstractManifold, + admo::AbstractDecoratedManifoldObjective, + m, + X, + n, +) + return linearized_forward_operator(M, N, get_objective(admo, false), m, X, n) +end + function linearized_forward_operator!( M::AbstractManifold, N::AbstractManifold, @@ -299,6 +332,17 @@ function linearized_forward_operator!( apdmo.linearized_forward_operator!!(M, Y, m, X) return Y end +function linearized_forward_operator!( + M::AbstractManifold, + N::AbstractManifold, + Y, + admo::AbstractDecoratedManifoldObjective, + m, + X, + n, +) + return linearized_forward_operator!(M, N, Y, get_objective(admo, false), m, X, n) +end @doc raw""" q = forward_operator(M::AbstractManifold, N::AbstractManifold, apdmo::AbstractPrimalDualManifoldObjective, p) @@ -338,6 +382,12 @@ function forward_operator( apdmo.Λ!!(M, q, p) return q end +function forward_operator( + M::AbstractManifold, N::AbstractManifold, admo::AbstractDecoratedManifoldObjective, p +) + return forward_operator(M, N, get_objective(admo, false), p) +end + function forward_operator!( M::AbstractManifold, N::AbstractManifold, @@ -358,6 +408,11 @@ function forward_operator!( apdmo.Λ!!(M, q, p) return q end +function forward_operator!( + M::AbstractManifold, N::AbstractManifold, q, admo::AbstractDecoratedManifoldObjective, p +) + return forward_operator!(M, N, q, get_objective(admo, false), p) +end @doc raw""" X = adjoint_linearized_operator(N::AbstractManifold, apdmo::AbstractPrimalDualManifoldObjective, m, n, Y) @@ -382,7 +437,6 @@ function adjoint_linearized_operator!(tmp::TwoManifoldProblem, X, m, n, Y) get_manifold(tmp, 1), get_manifold(tmp, 2), X, get_objective(tmp), m, n, Y ) end - function adjoint_linearized_operator( ::AbstractManifold, N::AbstractManifold, @@ -405,6 +459,17 @@ function adjoint_linearized_operator( apdmo.adjoint_linearized_operator!!(N, X, m, n, Y) return X end +function adjoint_linearized_operator( + M::AbstractManifold, + N::AbstractManifold, + admo::AbstractDecoratedManifoldObjective, + m, + n, + Y, +) + return adjoint_linearized_operator(M, N, get_objective(admo, false), m, n, Y) +end + function adjoint_linearized_operator!( M::AbstractManifold, N::AbstractManifold, @@ -429,6 +494,17 @@ function adjoint_linearized_operator!( apdmo.adjoint_linearized_operator!!(N, X, m, n, Y) return X end +function adjoint_linearized_operator!( + M::AbstractManifold, + N::AbstractManifold, + X, + admo::AbstractDecoratedManifoldObjective, + m, + n, + Y, +) + return adjoint_linearized_operator!(M, N, X, get_objective(admo, false), m, n, Y) +end @doc raw""" AbstractPrimalDualSolverState diff --git a/src/plans/proximal_plan.jl b/src/plans/proximal_plan.jl index ebef6c2248..2bdfd1e3f5 100644 --- a/src/plans/proximal_plan.jl +++ b/src/plans/proximal_plan.jl @@ -86,6 +86,12 @@ function get_proximal_map( check_prox_number(length(mpo.proximal_maps!!), i) return mpo.proximal_maps!![i](M, λ, p) end +function get_proximal_map( + M::AbstractManifold, admo::AbstractDecoratedManifoldObjective, λ, p, i +) + return get_proximal_map(M, get_objective(admo, false), λ, p, i) +end + function get_proximal_map!( M::AbstractManifold, q, mpo::ManifoldProximalMapObjective{AllocatingEvaluation}, λ, p, i ) @@ -93,6 +99,11 @@ function get_proximal_map!( copyto!(M, q, mpo.proximal_maps!![i](M, λ, p)) return q end +function get_proximal_map!( + M::AbstractManifold, q, admo::AbstractDecoratedManifoldObjective, λ, p, i +) + return get_proximal_map!(M, q, get_objective(admo, false), λ, p, i) +end function get_proximal_map( M::AbstractManifold, mpo::ManifoldProximalMapObjective{InplaceEvaluation}, λ, p, i ) diff --git a/src/plans/solver_state.jl b/src/plans/solver_state.jl index 3b549b82bc..0b87977795 100644 --- a/src/plans/solver_state.jl +++ b/src/plans/solver_state.jl @@ -83,16 +83,23 @@ show(io::IO, rst::ReturnSolverState) = print(io, "ReturnSolverState($(rst.state) dispatch_state_decorator(::ReturnSolverState) = Val(true) """ - get_solver_return(O::AbstractManoptSolverState) + get_solver_return(s::AbstractManoptSolverState) + get_solver_return(o::AbstractManifoldObjective, s::AbstractManoptSolverState) -determine the result value of a call to a solver. By default this returns the same as [`get_solver_result`](@ref), +determine the result value of a call to a solver. +By default this returns the same as [`get_solver_result`](@ref), i.e. the last iterate or (approximate) minimizer. - get_solver_return(O::ReturnSolverState) + get_solver_return(s::ReturnSolverState) + get_solver_return(o::AbstractManifoldObjective, s::ReturnSolverState) return the internally stored state of the [`ReturnSolverState`](@ref) instead of the minimizer. This means that when the state are decorated like this, the user still has to call [`get_solver_result`](@ref) on the internal state separately. + + get_solver_return(o::ReturnManifoldObjective, s::AbstractManoptSolverState) + +return both the objective and the state as a tuple. """ function get_solver_return(s::AbstractManoptSolverState) return _get_solver_return(s, dispatch_state_decorator(s)) @@ -100,6 +107,18 @@ end _get_solver_return(s::AbstractManoptSolverState, ::Val{false}) = get_solver_result(s) _get_solver_return(s::AbstractManoptSolverState, ::Val{true}) = get_solver_return(s.state) get_solver_return(s::ReturnSolverState) = s.state +function get_solver_return(o::AbstractManifoldObjective, s::AbstractManoptSolverState) + #resolve objevctive first + return _get_solver_return(o, s, dispatch_objective_decorator(o)) +end +#carefully undecorate both and check whether a solver/obejctive return happens +function _get_solver_return(o::AbstractManifoldObjective, s, ::Val{true}) + return get_solver_return(get_objective(o, false), s) +end +_get_solver_return(::AbstractManifoldObjective, s, ::Val{false}) = get_solver_return(s) +function get_solver_return(o::ReturnManifoldObjective, s::AbstractManoptSolverState) + return o.objective, get_solver_return(s) +end @doc raw""" get_state(s::AbstractManoptSolverState) @@ -174,15 +193,35 @@ _set_iterate!(s::AbstractManoptSolverState, M, p, ::Val{true}) = set_iterate!(s. """ get_solver_result(ams::AbstractManoptSolverState) + get_solver_result(tos::Tuple{AbstractManifoldObjective,AbstractManoptSolverState}) + get_solver_result(o::AbstractManifoldObjective, s::AbstractManoptSolverState) Return the final result after all iterations that is stored within the [`AbstractManoptSolverState`](@ref) `ams`, which was modified during the iterations. + +For the case the objective is passed as well, but default, the objective is ignored, +and the solver result for the state is called. """ function get_solver_result(s::AbstractManoptSolverState) - return get_solver_result(s, dispatch_state_decorator(s)) + return _get_solver_result(s, dispatch_state_decorator(s)) +end +function get_solver_result( + tos::Tuple{<:AbstractManifoldObjective,<:AbstractManoptSolverState} +) + return get_solver_result(tos...) end -get_solver_result(s::AbstractManoptSolverState, ::Val{false}) = get_iterate(s) -get_solver_result(s::AbstractManoptSolverState, ::Val{true}) = get_solver_result(s.state) +function get_solver_result(tos::Tuple{<:AbstractManifoldObjective,S}) where {S} + return tos[2] +end +function get_solver_result(::AbstractManifoldObjective, s::AbstractManoptSolverState) + return get_solver_result(s) +end +# if the second one is anything else we assume it is a point/result -> return that +function get_solver_result(::AbstractManifoldObjective, s) + return s +end +_get_solver_result(s::AbstractManoptSolverState, ::Val{false}) = get_iterate(s) +_get_solver_result(s::AbstractManoptSolverState, ::Val{true}) = get_solver_result(s.state) """ struct PointStorageKey{key} end @@ -235,7 +274,7 @@ _storage_copy_point(::AbstractManifold, p::Number) = StorageRef(p) Make a copy of tangent vector `X` from manifold `M` for storage in [`StoreStateAction`](@ref). """ -_storage_copy_vector(M::AbstractManifold, X) = copy(M, SA_F64[], X) +_storage_copy_vector(M::AbstractManifold, X) = copy(M, X) _storage_copy_vector(::AbstractManifold, X::Number) = StorageRef(X) @doc raw""" @@ -392,11 +431,9 @@ end function (a::StoreStateAction)( amp::AbstractManoptProblem, s::AbstractManoptSolverState, i::Int ) - #update values (maybe only once) - if !a.once || a.last_stored != i - update_storage!(a, amp, s) - end - return a.last_stored = i + (!a.once || a.last_stored != i) && (update_storage!(a, amp, s)) + a.last_stored = i + return a end """ @@ -564,3 +601,14 @@ end function get_count(ams::AbstractManoptSolverState, v::Val{:Iterations}) return get_count(ams.stop, v) end + +# in general, ignore printing the objective by default +function show(io::IO, t::Tuple{<:AbstractManifoldObjective,<:AbstractManoptSolverState}) + return print(io, "$(t[2])") +end +# for decorated ons, default: pass down +function show( + io::IO, t::Tuple{<:AbstractDecoratedManifoldObjective,<:AbstractManoptSolverState} +) + return show(io, (get_objective(t[1], false), t[2])) +end diff --git a/src/plans/stepsize.jl b/src/plans/stepsize.jl index 30ff265231..90604d7994 100644 --- a/src/plans/stepsize.jl +++ b/src/plans/stepsize.jl @@ -278,9 +278,9 @@ function show(io::IO, als::ArmijoLinesearch) io, """ ArmijoLineseach() with keyword parameters - * initial_stepsize = $(als.initial_stepsize) - * retraction_method = $(als.retraction_method) - * contraction_factor = $(als.contraction_factor) + * initial_stepsize = $(als.initial_stepsize) + * retraction_method = $(als.retraction_method) + * contraction_factor = $(als.contraction_factor) * sufficient_decrease = $(als.sufficient_decrease)""", ) end diff --git a/src/plans/stochastic_gradient_plan.jl b/src/plans/stochastic_gradient_plan.jl index f5344070d9..e09052454c 100644 --- a/src/plans/stochastic_gradient_plan.jl +++ b/src/plans/stochastic_gradient_plan.jl @@ -84,6 +84,10 @@ function get_gradients( ) where {TC} return [grad_i(M, p) for grad_i in sgo.gradient!!] end +function get_gradients(M::AbstractManifold, admo::AbstractDecoratedManifoldObjective, p) + return get_gradients(M, get_objective(admo, false), p) +end + function get_gradients!( M::AbstractManifold, X, @@ -102,6 +106,10 @@ function get_gradients!( copyto!.(Ref(M), X, [grad_i(M, p) for grad_i in sgo.gradient!!]) return X end +function get_gradients!(M::AbstractManifold, X, admo::AbstractDecoratedManifoldObjective, p) + return get_gradients!(M, X, get_objective(admo, false), p) +end + function get_gradients( ::AbstractManifold, ::ManifoldStochasticGradientObjective{InplaceEvaluation,TC,<:Function}, @@ -182,6 +190,10 @@ function get_gradient( X = zero_vector(M, p) return get_gradient!(M, X, sgo, p, k) end +function get_gradient(M::AbstractManifold, admo::AbstractDecoratedManifoldObjective, p, k) + return get_gradient(M, get_objective(admo, false), p, k) +end + function get_gradient!( M::AbstractManifold, X, @@ -222,6 +234,12 @@ function get_gradient!( ) where {TC} return sgo.gradient!![k](M, X, p) end +function get_gradient!( + M::AbstractManifold, X, admo::AbstractDecoratedManifoldObjective, p, k +) + return get_gradient!(M, X, get_objective(admo, false), p, k) +end + # Passdown from problem function get_gradient(mp::AbstractManoptProblem, p, k) return get_gradient(get_manifold(mp), get_objective(mp), p, k) diff --git a/src/plans/subgradient_plan.jl b/src/plans/subgradient_plan.jl index b98ff1ae27..f66daaf598 100644 --- a/src/plans/subgradient_plan.jl +++ b/src/plans/subgradient_plan.jl @@ -63,6 +63,10 @@ function get_subgradient( X = zero_vector(M, p) return sgo.subgradient!!(M, X, p) end +function get_subgradient(M::AbstractManifold, admo::AbstractDecoratedManifoldObjective, p) + return get_subgradient(M, get_objective(admo, false), p) +end + function get_subgradient!( M::AbstractManifold, X, sgo::ManifoldSubgradientObjective{AllocatingEvaluation}, p ) @@ -75,3 +79,8 @@ function get_subgradient!( sgo.subgradient!!(M, X, p) return X end +function get_subgradient!( + M::AbstractManifold, X, admo::AbstractDecoratedManifoldObjective, p +) + return get_subgradient!(M, X, get_objective(admo, false), p) +end diff --git a/src/solvers/ChambollePock.jl b/src/solvers/ChambollePock.jl index 3e37c23274..b1e4eca3d4 100644 --- a/src/solvers/ChambollePock.jl +++ b/src/solvers/ChambollePock.jl @@ -36,12 +36,15 @@ If you activate these to be different from the default identity, you have to pro `p.Λ` for the algorithm to work (which might be `missing` in the linearized case). # Constructor - ChambollePockState(M::AbstractManifold, + + ChambollePockState(M::AbstractManifold, N::AbstractManifold, m::P, n::Q, p::P, X::T, primal_stepsize::Float64, dual_stepsize::Float64; kwargs... ) -where all other fields from above are keyword arguments with their default values given in brackets, -as well as `N=TangentBundle(M)` + +where all other fields from above are keyword arguments with their default values given in brackets. + +if `Manifolds.jl` is loaded, `N` is also a keyword argument and set to `TangentBundle(M)` by default. """ mutable struct ChambollePockState{ P, @@ -73,65 +76,62 @@ mutable struct ChambollePockState{ inverse_retraction_method_dual::IRM_Dual vector_transport_method::VTM vector_transport_method_dual::VTM_Dual - - function ChambollePockState( - M::AbstractManifold, - m::P, - n::Q, - p::P, - X::T; - N=TangentBundle(M), - primal_stepsize::Float64=1 / sqrt(8), - dual_stepsize::Float64=1 / sqrt(8), - acceleration::Float64=0.0, - relaxation::Float64=1.0, - relax::Symbol=:primal, - stopping_criterion::StoppingCriterion=StopAfterIteration(300), - variant::Symbol=:exact, - update_primal_base::Union{Function,Missing}=missing, - update_dual_base::Union{Function,Missing}=missing, - retraction_method::RM=default_retraction_method(M, typeof(p)), - inverse_retraction_method::IRM=default_inverse_retraction_method(M, typeof(p)), - inverse_retraction_method_dual::IRM_Dual=default_inverse_retraction_method( - N, typeof(p) - ), - vector_transport_method::VTM=default_vector_transport_method(M, typeof(n)), - vector_transport_method_dual::VTM_Dual=default_vector_transport_method( - N, typeof(n) - ), - ) where { - P, - Q, - T, - RM<:AbstractRetractionMethod, - IRM<:AbstractInverseRetractionMethod, - IRM_Dual<:AbstractInverseRetractionMethod, - VTM<:AbstractVectorTransportMethod, - VTM_Dual<:AbstractVectorTransportMethod, - } - return new{P,Q,T,RM,IRM,IRM_Dual,VTM,VTM_Dual}( - m, - n, - p, - copy(M, p), - X, - copy(N, X), - primal_stepsize, - dual_stepsize, - acceleration, - relaxation, - relax, - stopping_criterion, - variant, - update_primal_base, - update_dual_base, - retraction_method, - inverse_retraction_method, - inverse_retraction_method_dual, - vector_transport_method, - vector_transport_method_dual, - ) - end +end +function Manopt.ChambollePockState( + M::AbstractManifold, + N::AbstractManifold, + m::P, + n::Q, + p::P, + X::T; + primal_stepsize::Float64=1 / sqrt(8), + dual_stepsize::Float64=1 / sqrt(8), + acceleration::Float64=0.0, + relaxation::Float64=1.0, + relax::Symbol=:primal, + stopping_criterion::StoppingCriterion=StopAfterIteration(300), + variant::Symbol=:exact, + update_primal_base::Union{Function,Missing}=missing, + update_dual_base::Union{Function,Missing}=missing, + retraction_method::RM=default_retraction_method(M, typeof(p)), + inverse_retraction_method::IRM=default_inverse_retraction_method(M, typeof(p)), + inverse_retraction_method_dual::IRM_Dual=default_inverse_retraction_method( + N, typeof(p) + ), + vector_transport_method::VTM=default_vector_transport_method(M, typeof(n)), + vector_transport_method_dual::VTM_Dual=default_vector_transport_method(N, typeof(n)), +) where { + P, + Q, + T, + RM<:AbstractRetractionMethod, + IRM<:AbstractInverseRetractionMethod, + IRM_Dual<:AbstractInverseRetractionMethod, + VTM<:AbstractVectorTransportMethod, + VTM_Dual<:AbstractVectorTransportMethod, +} + return ChambollePockState{P,Q,T,RM,IRM,IRM_Dual,VTM,VTM_Dual}( + m, + n, + p, + copy(M, p), + X, + copy(N, X), + primal_stepsize, + dual_stepsize, + acceleration, + relaxation, + relax, + stopping_criterion, + variant, + update_primal_base, + update_dual_base, + retraction_method, + inverse_retraction_method, + inverse_retraction_method_dual, + vector_transport_method, + vector_transport_method_dual, + ) end function show(io::IO, cps::ChambollePockState) i = get_count(cps, :Iterations) @@ -315,7 +315,7 @@ function ChambollePock!( ) dpdmo = decorate_objective!(M, pdmo; kwargs...) tmp = TwoManifoldProblem(M, N, dpdmo) - o = ChambollePockState( + cps = ChambollePockState( M, m, n, @@ -335,8 +335,9 @@ function ChambollePock!( inverse_retraction_method=inverse_retraction_method, vector_transport_method=vector_transport_method, ) - o = decorate_state!(o; kwargs...) - return get_solver_return(solve!(tmp, o)) + dcps = decorate_state!(cps; kwargs...) + solve!(tmp, dcps) + return get_solver_return(get_objective(tmp), dcps) end function initialize_solver!(::TwoManifoldProblem, ::ChambollePockState) end diff --git a/src/solvers/DouglasRachford.jl b/src/solvers/DouglasRachford.jl index 3aef0cd8da..5b5967ec45 100644 --- a/src/solvers/DouglasRachford.jl +++ b/src/solvers/DouglasRachford.jl @@ -171,8 +171,8 @@ function DouglasRachford( return (typeof(q) == typeof(rs)) ? rs[] : rs end function DouglasRachford( - M::AbstractManifold, mpo::ManifoldProximalMapObjective, p; kwargs... -) + M::AbstractManifold, mpo::O, p; kwargs... +) where {O<:Union{ManifoldProximalMapObjective,AbstractDecoratedManifoldObjective}} q = copy(M, p) return DouglasRachford!(M, mpo, q; kwargs...) end @@ -227,7 +227,7 @@ function DouglasRachford!( end function DouglasRachford!( M::AbstractManifold, - mpo::ManifoldProximalMapObjective, + mpo::O, p; λ::Tλ=(iter) -> 1.0, α::Tα=(iter) -> 0.9, @@ -237,14 +237,15 @@ function DouglasRachford!( StopAfterIteration(200), StopWhenChangeLess(10.0^-5) ), kwargs..., #especially may contain decorator options -) where {Tλ,Tα,TR} +) where {Tλ,Tα,TR,O<:Union{ManifoldProximalMapObjective,AbstractDecoratedManifoldObjective}} dmpo = decorate_objective!(M, mpo; kwargs...) dmp = DefaultManoptProblem(M, dmpo) drs = DouglasRachfordState( M, p; λ=λ, α=α, R=R, stopping_criterion=stopping_criterion, parallel=parallel > 0 ) ddrs = decorate_state!(drs; kwargs...) - return get_solver_return(solve!(dmp, ddrs)) + solve!(dmp, ddrs) + return get_solver_return(get_objective(dmp), ddrs) end # # An internal function that turns more than 2 proxes into a parallel variant diff --git a/src/solvers/FrankWolfe.jl b/src/solvers/FrankWolfe.jl index bf828c2079..47d875ccae 100644 --- a/src/solvers/FrankWolfe.jl +++ b/src/solvers/FrankWolfe.jl @@ -201,8 +201,8 @@ function Frank_Wolfe_method( return (typeof(q) == typeof(rs)) ? rs[] : rs end function Frank_Wolfe_method( - M::AbstractManifold, mgo::ManifoldGradientObjective, p; kwargs... -) + M::AbstractManifold, mgo::O, p; kwargs... +) where {O<:Union{ManifoldGradientObjective,AbstractDecoratedManifoldObjective}} q = copy(M, p) return Frank_Wolfe_method!(M, mgo, q; kwargs...) end @@ -229,7 +229,7 @@ function Frank_Wolfe_method!( end function Frank_Wolfe_method!( M::AbstractManifold, - mgo::AbstractManifoldGradientObjective, + mgo::O, p; initial_vector=zero_vector(M, p), evaluation=AllocatingEvaluation(), @@ -261,7 +261,11 @@ function Frank_Wolfe_method!( ) end, kwargs..., #collect rest -) where {TStop<:StoppingCriterion,TStep<:Stepsize} +) where { + TStop<:StoppingCriterion, + TStep<:Stepsize, + O<:Union{ManifoldGradientObjective,AbstractDecoratedManifoldObjective}, +} dmgo = decorate_objective!(M, mgo; kwargs...) dmp = DefaultManoptProblem(M, dmgo) fws = FrankWolfeState( @@ -274,8 +278,9 @@ function Frank_Wolfe_method!( stepsize=stepsize, stopping_criterion=stopping_criterion, ) - fws = decorate_state!(fws; kwargs...) - return get_solver_return(solve!(dmp, fws)) + dfws = decorate_state!(fws; kwargs...) + solve!(dmp, dfws) + return get_solver_return(get_objective(dmp), dfws) end function initialize_solver!(amp::AbstractManoptProblem, fws::FrankWolfeState) get_gradient!(amp, fws.X, fws.p) diff --git a/src/solvers/LevenbergMarquardt.jl b/src/solvers/LevenbergMarquardt.jl index 283729e797..54167ac927 100644 --- a/src/solvers/LevenbergMarquardt.jl +++ b/src/solvers/LevenbergMarquardt.jl @@ -108,8 +108,8 @@ function LevenbergMarquardt( return LevenbergMarquardt(M, nlso, p; evaluation=evaluation, kwargs...) end function LevenbergMarquardt( - M::AbstractManifold, nlso::NonlinearLeastSquaresObjective, p; kwargs... -) + M::AbstractManifold, nlso::O, p; kwargs... +) where {O<:Union{NonlinearLeastSquaresObjective,AbstractDecoratedManifoldObjective}} q = copy(M, p) return LevenbergMarquardt!(M, nlso, q; kwargs...) end @@ -159,7 +159,7 @@ function LevenbergMarquardt!( end function LevenbergMarquardt!( M::AbstractManifold, - nlso::NonlinearLeastSquaresObjective, + nlso::O, p; retraction_method::AbstractRetractionMethod=default_retraction_method(M, typeof(p)), stopping_criterion::StoppingCriterion=StopAfterIteration(200) | @@ -168,7 +168,7 @@ function LevenbergMarquardt!( debug=[DebugWarnIfCostIncreases()], expect_zero_residual::Bool=false, kwargs..., #collect rest -) +) where {O<:Union{NonlinearLeastSquaresObjective,AbstractDecoratedManifoldObjective}} i_nlso = get_objective(nlso) # undeecorate – for safety dnlso = decorate_objective!(M, nlso; kwargs...) nlsp = DefaultManoptProblem(M, dnlso) @@ -181,8 +181,9 @@ function LevenbergMarquardt!( retraction_method=retraction_method, expect_zero_residual=expect_zero_residual, ) - lms = decorate_state!(lms; debug=debug, kwargs...) - return get_solver_return(solve!(nlsp, lms)) + dlms = decorate_state!(lms; debug=debug, kwargs...) + solve!(nlsp, dlms) + return get_solver_return(get_objective(nlsp), dlms) end # # Solver functions diff --git a/src/solvers/NelderMead.jl b/src/solvers/NelderMead.jl index 3a1527862d..2c2a1592ee 100644 --- a/src/solvers/NelderMead.jl +++ b/src/solvers/NelderMead.jl @@ -216,11 +216,8 @@ function NelderMead(M::AbstractManifold, f, population::NelderMeadSimplex; kwarg return NelderMead(M, mco, population; kwargs...) end function NelderMead( - M::AbstractManifold, - mco::AbstractManifoldCostObjective, - population::NelderMeadSimplex; - kwargs..., -) + M::AbstractManifold, mco::O, population::NelderMeadSimplex; kwargs... +) where {O<:Union{AbstractManifoldCostObjective,AbstractDecoratedManifoldObjective}} res_population = NelderMeadSimplex(copy.(Ref(M), population.pts)) return NelderMead!(M, mco, res_population; kwargs...) end @@ -248,7 +245,7 @@ function NelderMead!(M::AbstractManifold, f, population::NelderMeadSimplex; kwar end function NelderMead!( M::AbstractManifold, - mco::AbstractManifoldCostObjective, + mco::O, population::NelderMeadSimplex; stopping_criterion::StoppingCriterion=StopAfterIteration(2000) | StopWhenPopulationConcentrated(), @@ -263,7 +260,7 @@ function NelderMead!( M, eltype(population.pts) ), kwargs..., #collect rest -) +) where {O<:Union{AbstractManifoldCostObjective,AbstractDecoratedManifoldObjective}} dmco = decorate_objective!(M, mco; kwargs...) mp = DefaultManoptProblem(M, dmco) s = NelderMeadState( diff --git a/src/solvers/alternating_gradient_descent.jl b/src/solvers/alternating_gradient_descent.jl index 647adc8094..d212a84caa 100644 --- a/src/solvers/alternating_gradient_descent.jl +++ b/src/solvers/alternating_gradient_descent.jl @@ -146,6 +146,9 @@ function default_stepsize(M::AbstractManifold, ::Type{AlternatingGradientDescent return ArmijoLinesearch(M) end +function alternating_gradient_descent end +function alternating_gradient_descent! end + @doc raw""" alternating_gradient_descent(M::ProductManifold, f, grad_f, p=rand(M)) alternating_gradient_descent(M::ProductManifold, ago::ManifoldAlternatingGradientObjective, p) @@ -180,7 +183,7 @@ usually the obtained (approximate) minimizer, see [`get_solver_return`](@ref) fo !!! note - This Problem requires the `ProductManifold` from `Manifolds.jl`, so `Manifolds.jl` to be loaded. + This Problem requires the `ProductManifold` from `Manifolds.jl`, so `Manifolds.jl` needs to be loaded. !!! note @@ -189,24 +192,8 @@ usually the obtained (approximate) minimizer, see [`get_solver_return`](@ref) fo the `i`th components gradient is computed / returned. """ -alternating_gradient_descent(::ProductManifold, args...; kwargs...) -function alternating_gradient_descent( - M::ProductManifold, - f, - grad_f::Union{TgF,AbstractVector{<:TgF}}, - p=rand(M); - evaluation::AbstractEvaluationType=AllocatingEvaluation(), - kwargs..., -) where {TgF} - ago = ManifoldAlternatingGradientObjective(f, grad_f; evaluation=evaluation) - return alternating_gradient_descent(M, ago, p; evaluation=evaluation, kwargs...) -end -function alternating_gradient_descent( - M::ProductManifold, ago::ManifoldAlternatingGradientObjective, p; kwargs... -) - q = copy(M, p) - return alternating_gradient_descent!(M, ago, q; kwargs...) -end +alternating_gradient_descent(::AbstractManifold, args...; kwargs...) + @doc raw""" alternating_gradient_descent!(M::ProductManifold, f, grad_f, p) alternating_gradient_descent!(M::ProductManifold, ago::ManifoldAlternatingGradientObjective, p) @@ -225,47 +212,8 @@ you can also pass a [`ManifoldAlternatingGradientObjective`](@ref) `ago` contain for all optional parameters, see [`alternating_gradient_descent`](@ref). """ -alternating_gradient_descent!(M::ProductManifold, args...; kwargs...) - -function alternating_gradient_descent!( - M::ProductManifold, - f, - grad_f::Union{TgF,AbstractVector{<:TgF}}, - p; - evaluation::AbstractEvaluationType=AllocatingEvaluation(), - kwargs..., -) where {TgF} - agmo = ManifoldAlternatingGradientObjective(f, grad_f; evaluation=evaluation) - return alternating_gradient_descent!(M, agmo, p; evaluation=evaluation, kwargs...) -end -function alternating_gradient_descent!( - M::ProductManifold, - agmo::ManifoldAlternatingGradientObjective, - p; - inner_iterations::Int=5, - stopping_criterion::StoppingCriterion=StopAfterIteration(100) | - StopWhenGradientNormLess(1e-9), - stepsize::Stepsize=default_stepsize(M, AlternatingGradientDescentState), - order_type::Symbol=:Linear, - order=collect(1:length(M.manifolds)), - retraction_method::AbstractRetractionMethod=default_retraction_method(M, typeof(p)), - kwargs..., -) - dagmo = decorate_objective!(M, agmo; kwargs...) - dmp = DefaultManoptProblem(M, dagmo) - agds = AlternatingGradientDescentState( - M, - p; - inner_iterations=inner_iterations, - stopping_criterion=stopping_criterion, - stepsize=stepsize, - order_type=order_type, - order=order, - retraction_method=retraction_method, - ) - agds = decorate_state!(agds; kwargs...) - return get_solver_return(solve!(dmp, agds)) -end +alternating_gradient_descent!(M::AbstractManifold, args...; kwargs...) + function initialize_solver!( amp::AbstractManoptProblem, agds::AlternatingGradientDescentState ) diff --git a/src/solvers/augmented_Lagrangian_method.jl b/src/solvers/augmented_Lagrangian_method.jl index ba92b277e9..004029ffa7 100644 --- a/src/solvers/augmented_Lagrangian_method.jl +++ b/src/solvers/augmented_Lagrangian_method.jl @@ -261,8 +261,8 @@ function augmented_Lagrangian_method( return augmented_Lagrangian_method!(M, cmo, q; evaluation=evaluation, kwargs...) end function augmented_Lagrangian_method( - M::AbstractManifold, cmo::ConstrainedManifoldObjective, p=rand(M); kwargs... -) + M::AbstractManifold, cmo::O, p=rand(M); kwargs... +) where {O<:Union{ConstrainedManifoldObjective,AbstractDecoratedManifoldObjective}} q = copy(M, p) return augmented_Lagrangian_method!(M, cmo, q; kwargs...) end @@ -319,7 +319,7 @@ function augmented_Lagrangian_method!( end function augmented_Lagrangian_method!( M::AbstractManifold, - cmo::ConstrainedManifoldObjective, + cmo::O, p; evaluation=AllocatingEvaluation(), ϵ::Real=1e-3, @@ -360,7 +360,7 @@ function augmented_Lagrangian_method!( StopWhenSmallerOrEqual(:ϵ, ϵ_min) & StopWhenChangeLess(1e-10) ), kwargs..., -) +) where {O<:Union{ConstrainedManifoldObjective,AbstractDecoratedManifoldObjective}} alms = AugmentedLagrangianMethodState( M, cmo, @@ -383,7 +383,8 @@ function augmented_Lagrangian_method!( dcmo = decorate_objective!(M, cmo; kwargs...) mp = DefaultManoptProblem(M, dcmo) alms = decorate_state!(alms; kwargs...) - return get_solver_return(solve!(mp, alms)) + solve!(mp, alms) + return get_solver_return(get_objective(mp), alms) end # @@ -410,14 +411,14 @@ function step_solver!(mp::AbstractManoptProblem, alms::AugmentedLagrangianMethod # update multipliers cost_ineq = get_inequality_constraints(mp, alms.p) - n_ineq_constraint = size(cost_ineq, 1) + n_ineq_constraint = length(cost_ineq) alms.μ .= min.( ones(n_ineq_constraint) .* alms.μ_max, max.(alms.μ .+ alms.ρ .* cost_ineq, zeros(n_ineq_constraint)), ) cost_eq = get_equality_constraints(mp, alms.p) - n_eq_constraint = size(cost_eq, 1) + n_eq_constraint = length(cost_eq) alms.λ = min.( ones(n_eq_constraint) .* alms.λ_max, diff --git a/src/solvers/conjugate_gradient_descent.jl b/src/solvers/conjugate_gradient_descent.jl index f69efbfa19..427b81d2eb 100644 --- a/src/solvers/conjugate_gradient_descent.jl +++ b/src/solvers/conjugate_gradient_descent.jl @@ -111,8 +111,8 @@ function conjugate_gradient_descent( return (typeof(q) == typeof(rs)) ? rs[] : rs end function conjugate_gradient_descent( - M::AbstractManifold, mgo::ManifoldGradientObjective, p=rand(M); kwargs... -) + M::AbstractManifold, mgo::O, p=rand(M); kwargs... +) where {O<:Union{ManifoldGradientObjective,AbstractDecoratedManifoldObjective}} q = copy(M, p) return conjugate_gradient_descent!(M, mgo, q; kwargs...) end @@ -154,7 +154,7 @@ function conjugate_gradient_descent!( end function conjugate_gradient_descent!( M::AbstractManifold, - mgo::ManifoldGradientObjective, + mgo::O, p; coefficient::DirectionUpdateRule=ConjugateDescentCoefficient(), retraction_method::AbstractRetractionMethod=default_retraction_method(M, typeof(p)), @@ -166,7 +166,7 @@ function conjugate_gradient_descent!( ), vector_transport_method=default_vector_transport_method(M, typeof(p)), kwargs..., -) +) where {O<:Union{ManifoldGradientObjective,AbstractDecoratedManifoldObjective}} dmgo = decorate_objective!(M, mgo; kwargs...) dmp = DefaultManoptProblem(M, dmgo) X = zero_vector(M, p) @@ -180,8 +180,9 @@ function conjugate_gradient_descent!( vector_transport_method, X, ) - cgs = decorate_state!(cgs; kwargs...) - return get_solver_return(solve!(dmp, cgs)) + dcgs = decorate_state!(cgs; kwargs...) + solve!(dmp, dcgs) + return get_solver_return(get_objective(dmp), dcgs) end function initialize_solver!(amp::AbstractManoptProblem, cgs::ConjugateGradientDescentState) cgs.X = get_gradient(amp, cgs.p) diff --git a/src/solvers/cyclic_proximal_point.jl b/src/solvers/cyclic_proximal_point.jl index a08a095ac2..c75c28bb3f 100644 --- a/src/solvers/cyclic_proximal_point.jl +++ b/src/solvers/cyclic_proximal_point.jl @@ -81,8 +81,8 @@ function cyclic_proximal_point( return (typeof(q) == typeof(rs)) ? rs[] : rs end function cyclic_proximal_point( - M::AbstractManifold, mpo::ManifoldProximalMapObjective, p; kwargs... -) + M::AbstractManifold, mpo::O, p; kwargs... +) where {O<:Union{ManifoldProximalMapObjective,AbstractDecoratedManifoldObjective}} q = copy(M, p) return cyclic_proximal_point!(M, mpo, q; kwargs...) end @@ -118,7 +118,7 @@ function cyclic_proximal_point!( end function cyclic_proximal_point!( M::AbstractManifold, - mpo::ManifoldProximalMapObjective, + mpo::O, p; evaluation_order::Symbol=:Linear, stopping_criterion::StoppingCriterion=StopWhenAny( @@ -126,14 +126,15 @@ function cyclic_proximal_point!( ), λ=i -> 1 / i, kwargs..., -) +) where {O<:Union{ManifoldProximalMapObjective,AbstractDecoratedManifoldObjective}} dmpo = decorate_objective!(M, mpo; kwargs...) dmp = DefaultManoptProblem(M, dmpo) cpps = CyclicProximalPointState( M, p; stopping_criterion=stopping_criterion, λ=λ, evaluation_order=evaluation_order ) - cpps = decorate_state!(cpps; kwargs...) - return get_solver_return(solve!(dmp, cpps)) + dcpps = decorate_state!(cpps; kwargs...) + solve!(dmp, dcpps) + return get_solver_return(get_objective(dmp), dcpps) end function initialize_solver!(amp::AbstractManoptProblem, cpps::CyclicProximalPointState) c = length(get_objective(amp).proximal_maps!!) diff --git a/src/solvers/difference-of-convex-proximal-point.jl b/src/solvers/difference-of-convex-proximal-point.jl index 12335c50a6..ae2ea127c8 100644 --- a/src/solvers/difference-of-convex-proximal-point.jl +++ b/src/solvers/difference-of-convex-proximal-point.jl @@ -269,8 +269,10 @@ function difference_of_convex_proximal_point( end function difference_of_convex_proximal_point( - M::AbstractManifold, mdcpo::ManifoldDifferenceOfConvexProximalObjective, p; kwargs... -) + M::AbstractManifold, mdcpo::O, p; kwargs... +) where { + O<:Union{ManifoldDifferenceOfConvexProximalObjective,AbstractDecoratedManifoldObjective} +} q = copy(M, p) return difference_of_convex_proximal_point!(M, mdcpo, q; kwargs...) end @@ -310,7 +312,7 @@ function difference_of_convex_proximal_point!( end function difference_of_convex_proximal_point!( M::AbstractManifold, - mdcpo::ManifoldDifferenceOfConvexProximalObjective, + mdcpo::O, p; g=nothing, grad_g=nothing, @@ -360,7 +362,9 @@ function difference_of_convex_proximal_point!( ) end, kwargs..., -) +) where { + O<:Union{ManifoldDifferenceOfConvexProximalObjective,AbstractDecoratedManifoldObjective} +} # Check whether either the right defaults were provided or a sub_problen. if isnothing(sub_problem) error( @@ -389,7 +393,8 @@ function difference_of_convex_proximal_point!( λ=λ, ) ddcps = decorate_state!(dcps; kwargs...) - return get_solver_return(solve!(dmp, ddcps)) + solve!(dmp, ddcps) + return get_solver_return(get_objective(dmp), ddcps) end function initialize_solver!(::AbstractManoptProblem, dcps::DifferenceOfConvexProximalState) diff --git a/src/solvers/difference_of_convex_algorithm.jl b/src/solvers/difference_of_convex_algorithm.jl index 4f9a0cb12f..d4adb75e58 100644 --- a/src/solvers/difference_of_convex_algorithm.jl +++ b/src/solvers/difference_of_convex_algorithm.jl @@ -234,8 +234,8 @@ function difference_of_convex_algorithm( return (typeof(q) == typeof(rs)) ? rs[] : rs end function difference_of_convex_algorithm( - M::AbstractManifold, mdco::ManifoldDifferenceOfConvexObjective, p; kwargs... -) + M::AbstractManifold, mdco::O, p; kwargs... +) where {O<:Union{ManifoldDifferenceOfConvexObjective,AbstractDecoratedManifoldObjective}} q = copy(M, p) return difference_of_convex_algorithm!(M, mdco, p; kwargs...) end @@ -268,7 +268,7 @@ function difference_of_convex_algorithm!( end function difference_of_convex_algorithm!( M::AbstractManifold, - mdco::ManifoldDifferenceOfConvexObjective, + mdco::O, p; evaluation::AbstractEvaluationType=AllocatingEvaluation(), g=nothing, @@ -316,7 +316,7 @@ function difference_of_convex_algorithm!( ) end, kwargs..., #collect rest -) +) where {O<:Union{ManifoldDifferenceOfConvexObjective,AbstractDecoratedManifoldObjective}} dmdco = decorate_objective!(M, mdco; kwargs...) dmp = DefaultManoptProblem(M, dmdco) # For now only subsolvers - TODO closed form solution init here @@ -350,7 +350,8 @@ function difference_of_convex_algorithm!( ) end ddcs = decorate_state!(dcs; kwargs...) - return get_solver_return(solve!(dmp, ddcs)) + solve!(dmp, ddcs) + return get_solver_return(get_objective(dmp), ddcs) end function initialize_solver!(::AbstractManoptProblem, dcs::DifferenceOfConvexState) return dcs diff --git a/src/solvers/exact_penalty_method.jl b/src/solvers/exact_penalty_method.jl index 01e895245a..75b2dcf322 100644 --- a/src/solvers/exact_penalty_method.jl +++ b/src/solvers/exact_penalty_method.jl @@ -260,8 +260,8 @@ function exact_penalty_method( return (typeof(q) == typeof(rs)) ? rs[] : rs end function exact_penalty_method( - M::AbstractManifold, cmo::ConstrainedManifoldObjective, p=rand(M); kwargs... -) + M::AbstractManifold, cmo::O, p=rand(M); kwargs... +) where {O<:Union{ConstrainedManifoldObjective,AbstractDecoratedManifoldObjective}} q = copy(M, p) return exact_penalty_method!(M, cmo, q; kwargs...) end @@ -294,7 +294,7 @@ function exact_penalty_method!( end function exact_penalty_method!( M::AbstractManifold, - cmo::ConstrainedManifoldObjective, + cmo::O, p; evaluation=AllocatingEvaluation(), ϵ::Real=1e-3, @@ -334,7 +334,7 @@ function exact_penalty_method!( StopWhenSmallerOrEqual(:ϵ, ϵ_min) & StopWhenChangeLess(1e-10) ), kwargs..., -) +) where {O<:Union{ConstrainedManifoldObjective,AbstractDecoratedManifoldObjective}} emps = ExactPenaltyMethodState( M, p, @@ -353,7 +353,8 @@ function exact_penalty_method!( deco_o = decorate_objective!(M, cmo; kwargs...) dmp = DefaultManoptProblem(M, deco_o) epms = decorate_state!(emps; kwargs...) - return get_solver_return(solve!(dmp, epms)) + solve!(dmp, epms) + return get_solver_return(get_objective(dmp), epms) end # # Solver functions diff --git a/src/solvers/gradient_descent.jl b/src/solvers/gradient_descent.jl index 76b6076b19..8eccbc510c 100644 --- a/src/solvers/gradient_descent.jl +++ b/src/solvers/gradient_descent.jl @@ -146,8 +146,8 @@ the [`AbstractManifoldGradientObjective`](@ref) `gradient_objective` directly. If you provide the [`ManifoldGradientObjective`](@ref) directly, `evaluation` is ignored. -All other keyword arguments are passed to [`decorate_state!`](@ref) for decorators or -[`decorate_objective!`](@ref), respectively. +All other keyword arguments are passed to [`decorate_state!`](@ref) for state decorators or +[`decorate_objective!`](@ref) for objective, respectively. If you provide the [`ManifoldGradientObjective`](@ref) directly, these decorations can still be specified # Output @@ -186,7 +186,9 @@ function gradient_descent( #return just a number if the return type is the same as the type of q return (typeof(q) == typeof(rs)) ? rs[] : rs end -function gradient_descent(M::AbstractManifold, mgo::ManifoldGradientObjective, p; kwargs...) +function gradient_descent( + M::AbstractManifold, mgo::O, p; kwargs... +) where {O<:Union{AbstractManifoldGradientObjective,AbstractDecoratedManifoldObjective}} q = copy(M, p) return gradient_descent!(M, mgo, q; kwargs...) end @@ -228,7 +230,7 @@ function gradient_descent!( end function gradient_descent!( M::AbstractManifold, - mgo::AbstractManifoldGradientObjective, + mgo::O, p; retraction_method::AbstractRetractionMethod=default_retraction_method(M, typeof(p)), stepsize::Stepsize=default_stepsize( @@ -239,9 +241,9 @@ function gradient_descent!( debug=stepsize isa ConstantStepsize ? [DebugWarnIfCostIncreases()] : [], direction=IdentityUpdateRule(), kwargs..., #collect rest -) +) where {O<:Union{AbstractManifoldGradientObjective,AbstractDecoratedManifoldObjective}} dmgo = decorate_objective!(M, mgo; kwargs...) - mp = DefaultManoptProblem(M, dmgo) + dmp = DefaultManoptProblem(M, dmgo) s = GradientDescentState( M, p; @@ -250,8 +252,9 @@ function gradient_descent!( direction=direction, retraction_method=retraction_method, ) - s = decorate_state!(s; debug=debug, kwargs...) - return get_solver_return(solve!(mp, s)) + ds = decorate_state!(s; debug=debug, kwargs...) + solve!(dmp, ds) + return get_solver_return(get_objective(dmp), ds) end # # Solver functions diff --git a/src/solvers/particle_swarm.jl b/src/solvers/particle_swarm.jl index eecb9c74ae..90d115fa3b 100644 --- a/src/solvers/particle_swarm.jl +++ b/src/solvers/particle_swarm.jl @@ -251,11 +251,8 @@ function particle_swarm( end function particle_swarm( - M::AbstractManifold, - mco::AbstractManifoldCostObjective, - swarm::AbstractVector; - kwargs..., -) + M::AbstractManifold, mco::O, swarm::AbstractVector; kwargs... +) where {O<:Union{AbstractManifoldCostObjective,AbstractDecoratedManifoldObjective}} new_swarm = [copy(M, xi) for xi in swarm] return particle_swarm!(M, mco, new_swarm; kwargs...) end @@ -283,7 +280,7 @@ function particle_swarm!(M::AbstractManifold, f, swarm::AbstractVector; kwargs.. end function particle_swarm!( M::AbstractManifold, - mco::AbstractManifoldCostObjective, + mco::O, swarm::AbstractVector; velocity::AbstractVector=[rand(M; vector_at=y) for y in swarm], inertia::Real=0.65, @@ -299,7 +296,7 @@ function particle_swarm!( M, eltype(swarm) ), kwargs..., #collect rest -) +) where {O<:Union{AbstractManifoldCostObjective,AbstractDecoratedManifoldObjective}} dmco = decorate_objective!(M, mco; kwargs...) mp = DefaultManoptProblem(M, dmco) pss = ParticleSwarmState( @@ -314,8 +311,9 @@ function particle_swarm!( inverse_retraction_method=inverse_retraction_method, vector_transport_method=vector_transport_method, ) - o = decorate_state!(pss; kwargs...) - return get_solver_return(solve!(mp, o)) + dpss = decorate_state!(pss; kwargs...) + solve!(mp, dpss) + return get_solver_return(get_objective(mp), dpss) end # diff --git a/src/solvers/primal_dual_semismooth_Newton.jl b/src/solvers/primal_dual_semismooth_Newton.jl index da4450bc06..1d98cf32f0 100644 --- a/src/solvers/primal_dual_semismooth_Newton.jl +++ b/src/solvers/primal_dual_semismooth_Newton.jl @@ -154,8 +154,9 @@ function primal_dual_semismooth_Newton!( inverse_retraction_method=inverse_retraction_method, vector_transport_method=vector_transport_method, ) - pdsn = decorate_state!(pdsn; kwargs...) - return get_solver_return(solve!(tmp, pdsn)) + dpdsn = decorate_state!(pdsn; kwargs...) + solve!(tmp, dpdsn) + return get_solver_return(get_objective(tmp), dpdsn) end function initialize_solver!(::TwoManifoldProblem, ::PrimalDualSemismoothNewtonState) end diff --git a/src/solvers/quasi_Newton.jl b/src/solvers/quasi_Newton.jl index 8fe5094c9f..69c9e47625 100644 --- a/src/solvers/quasi_Newton.jl +++ b/src/solvers/quasi_Newton.jl @@ -220,7 +220,9 @@ function quasi_Newton( #return just a number if the return type is the same as the type of q return (typeof(q) == typeof(rs)) ? rs[] : rs end -function quasi_Newton(M::AbstractManifold, mgo::ManifoldGradientObjective, p; kwargs...) +function quasi_Newton( + M::AbstractManifold, mgo::O, p; kwargs... +) where {O<:Union{ManifoldGradientObjective,AbstractDecoratedManifoldObjective}} q = copy(M, p) return quasi_Newton!(M, mgo, q; kwargs...) end @@ -252,7 +254,7 @@ function quasi_Newton!( end function quasi_Newton!( M::AbstractManifold, - mgo::ManifoldGradientObjective, + mgo::O, p; cautious_update::Bool=false, cautious_function::Function=x -> x * 10^(-4), @@ -281,7 +283,7 @@ function quasi_Newton!( stopping_criterion::StoppingCriterion=StopAfterIteration(max(1000, memory_size)) | StopWhenGradientNormLess(1e-6), kwargs..., -) +) where {O<:Union{ManifoldGradientObjective,AbstractDecoratedManifoldObjective}} if memory_size >= 0 local_dir_upd = QuasiNewtonLimitedMemoryDirectionUpdate( M, @@ -319,8 +321,9 @@ function quasi_Newton!( retraction_method=retraction_method, vector_transport_method=vector_transport_method, ) - qns = decorate_state!(qns; kwargs...) - return get_solver_return(solve!(mp, qns)) + dqns = decorate_state!(qns; kwargs...) + solve!(mp, dqns) + return get_solver_return(get_objective(mp), dqns) end function initialize_solver!(p::AbstractManoptProblem, s::QuasiNewtonState) s.X = get_gradient(p, s.p) diff --git a/src/solvers/solver.jl b/src/solvers/solver.jl index bd85802503..10ec67ce58 100644 --- a/src/solvers/solver.jl +++ b/src/solvers/solver.jl @@ -64,9 +64,11 @@ decorate the [`AbstractManifoldObjective`](@ref)` o` with specific decorators. optional arguments provide necessary details on the decorators. A specific one is used to activate certain decorators. -* `cache` – (`missing`) currently only supports the [`SimpleCacheObjective`](@ref) - which is activated by either specifying the symbol `:Simple` or the tuple - (`:Simple, kwargs...`) to pass down keyword arguments +* `cache` – (`missing`) specify a cache. Currenlty `:Simple` is supported and `:LRU` if you + load `LRUCache.jl`. For this case a tuple specifying what to cache and how many can be provided, + i.e. `(:LRU, [:Cost, :Gradient], 10)`, where the number specifies the size of each cache. + and 10 is the default if one omits the last tuple entry +* `count` – (`missing`) specify calls to the objective to be called, see [`ManifoldCountObjective`](@ref) for the full list other keywords are ignored. @@ -75,9 +77,20 @@ other keywords are ignored. [`objective_cache_factory`](@ref) """ function decorate_objective!( - M::AbstractManifold, o::O; cache::Union{Missing,Symbol}=missing, kwargs... -) where {O<:AbstractManifoldObjective} - deco_o = ismissing(cache) ? o : objective_cache_factory(M, o, cache) + M::AbstractManifold, + o::O; + cache::Union{ + Missing,Symbol,Tuple{Symbol,<:AbstractArray},Tuple{Symbol,<:AbstractArray,P} + }=missing, + count::Union{Missing,AbstractVector{<:Symbol}}=missing, + return_objective=false, + kwargs..., +) where {O<:AbstractManifoldObjective,P} + # Order: First count _then_ cache, so that cache is the outer decorator and + # we only count _after_ cache misses + deco_o = ismissing(count) ? o : objective_count_factory(M, o, count) + deco_o = ismissing(cache) ? deco_o : objective_cache_factory(M, deco_o, cache) + deco_o = return_objective ? ReturnManifoldObjective(deco_o) : deco_o return deco_o end diff --git a/src/solvers/stochastic_gradient_descent.jl b/src/solvers/stochastic_gradient_descent.jl index f44344c9a6..d6dbca2515 100644 --- a/src/solvers/stochastic_gradient_descent.jl +++ b/src/solvers/stochastic_gradient_descent.jl @@ -202,8 +202,8 @@ function stochastic_gradient_descent( return (typeof(q) == typeof(rs)) ? rs[] : rs end function stochastic_gradient_descent( - M::AbstractManifold, msgo::ManifoldStochasticGradientObjective, p; kwargs... -) + M::AbstractManifold, msgo::O, p; kwargs... +) where {O<:Union{ManifoldStochasticGradientObjective,AbstractDecoratedManifoldObjective}} q = copy(M, p) return stochastic_gradient_descent!(M, msgo, q; kwargs...) end @@ -240,7 +240,7 @@ function stochastic_gradient_descent!( end function stochastic_gradient_descent!( M::AbstractManifold, - msgo::ManifoldStochasticGradientObjective, + msgo::O, p; direction::DirectionUpdateRule=StochasticGradient(zero_vector(M, p)), stopping_criterion::StoppingCriterion=StopAfterIteration(10000) | @@ -250,7 +250,7 @@ function stochastic_gradient_descent!( order_type::Symbol=:Random, retraction_method::AbstractRetractionMethod=default_retraction_method(M, typeof(p)), kwargs..., -) +) where {O<:Union{ManifoldStochasticGradientObjective,AbstractDecoratedManifoldObjective}} dmsgo = decorate_objective!(M, msgo; kwargs...) mp = DefaultManoptProblem(M, dmsgo) sgds = StochasticGradientDescentState( @@ -264,8 +264,9 @@ function stochastic_gradient_descent!( order=order, retraction_method=retraction_method, ) - sgds = decorate_state!(sgds; kwargs...) - return get_solver_return(solve!(mp, sgds)) + dsgds = decorate_state!(sgds; kwargs...) + solve!(mp, dsgds) + return get_solver_return(get_objective(mp), dsgds) end function initialize_solver!(::AbstractManoptProblem, s::StochasticGradientDescentState) s.k = 1 diff --git a/src/solvers/subgradient.jl b/src/solvers/subgradient.jl index 4932658d05..760fce04d9 100644 --- a/src/solvers/subgradient.jl +++ b/src/solvers/subgradient.jl @@ -149,8 +149,8 @@ function subgradient_method( return (typeof(q) == typeof(rs)) ? rs[] : rs end function subgradient_method( - M::AbstractManifold, sgo::ManifoldSubgradientObjective, p; kwargs... -) + M::AbstractManifold, sgo::O, p; kwargs... +) where {O<:Union{ManifoldSubgradientObjective,AbstractDecoratedManifoldObjective}} q = copy(M, p) return subgradient_method!(M, sgo, q; kwargs...) end @@ -189,13 +189,13 @@ function subgradient_method!( end function subgradient_method!( M::AbstractManifold, - sgo::ManifoldSubgradientObjective, + sgo::O, p; retraction_method::AbstractRetractionMethod=default_retraction_method(M, typeof(p)), stepsize::Stepsize=default_stepsize(M, SubGradientMethodState), stopping_criterion::StoppingCriterion=StopAfterIteration(5000), kwargs..., -) +) where {O<:Union{ManifoldSubgradientObjective,AbstractDecoratedManifoldObjective}} dsgo = decorate_objective!(M, sgo; kwargs...) mp = DefaultManoptProblem(M, dsgo) sgs = SubGradientMethodState( @@ -205,8 +205,9 @@ function subgradient_method!( stepsize=stepsize, retraction_method=retraction_method, ) - sgs = decorate_state!(sgs; kwargs...) - return get_solver_return(solve!(mp, sgs)) + dsgs = decorate_state!(sgs; kwargs...) + solve!(mp, dsgs) + return get_solver_return(get_objective(mp), dsgs) end function initialize_solver!(mp::AbstractManoptProblem, sgs::SubGradientMethodState) M = get_manifold(mp) diff --git a/src/solvers/truncated_conjugate_gradient_descent.jl b/src/solvers/truncated_conjugate_gradient_descent.jl index 5d8e140d26..706f676ada 100644 --- a/src/solvers/truncated_conjugate_gradient_descent.jl +++ b/src/solvers/truncated_conjugate_gradient_descent.jl @@ -518,8 +518,8 @@ end # Objetive -> Allocate and call ! # function truncated_conjugate_gradient_descent( - M::AbstractManifold, mho::ManifoldHessianObjective, p, X; kwargs... -) + M::AbstractManifold, mho::O, p, X; kwargs... +) where {O<:Union{ManifoldHessianObjective,AbstractDecoratedManifoldObjective}} q = copy(M, p) Y = copy(M, p, X) return truncated_conjugate_gradient_descent!(M, mho, q, Y; kwargs...) @@ -597,7 +597,7 @@ function truncated_conjugate_gradient_descent!( end function truncated_conjugate_gradient_descent!( M::AbstractManifold, - mho::ManifoldHessianObjective, + mho::O, p, X; trust_region_radius::Float64=injectivity_radius(M) / 4, @@ -613,7 +613,7 @@ function truncated_conjugate_gradient_descent!( StopWhenModelIncreased(), project!::Proj=copyto!, kwargs..., #collect rest -) where {Proj} +) where {Proj,O<:Union{ManifoldHessianObjective,AbstractDecoratedManifoldObjective}} dmho = decorate_objective!(M, mho; kwargs...) mp = DefaultManoptProblem(M, dmho) tcgs = TruncatedConjugateGradientState( @@ -627,8 +627,9 @@ function truncated_conjugate_gradient_descent!( stopping_criterion=stopping_criterion, (project!)=project!, ) - tcgs = decorate_state!(tcgs; kwargs...) - return get_solver_return(solve!(mp, tcgs)) + dtcgs = decorate_state!(tcgs; kwargs...) + solve!(mp, dtcgs) + return get_solver_return(get_objective(mp), dtcgs) end # # Deprecated - even kept old notation diff --git a/src/solvers/trust_regions.jl b/src/solvers/trust_regions.jl index dbafa490fd..ee1fcf4531 100644 --- a/src/solvers/trust_regions.jl +++ b/src/solvers/trust_regions.jl @@ -349,8 +349,8 @@ function trust_regions( ) end function trust_regions( - M::AbstractManifold, mho::ManifoldHessianObjective, p=rand(M); kwargs... -) + M::AbstractManifold, mho::O, p=rand(M); kwargs... +) where {O<:Union{ManifoldHessianObjective,AbstractDecoratedManifoldObjective}} q = copy(M, p) return trust_regions!(M, mho, q; kwargs...) end @@ -416,7 +416,7 @@ function trust_regions!( end function trust_regions!( M::AbstractManifold, - mho::ManifoldHessianObjective, + mho::O, p; retraction_method::AbstractRetractionMethod=default_retraction_method(M, typeof(p)), stopping_criterion::StoppingCriterion=StopAfterIteration(1000) | @@ -442,7 +442,7 @@ function trust_regions!( (project!)=project!, ), kwargs..., #collect rest -) where {Proj} +) where {Proj,O<:Union{ManifoldHessianObjective,AbstractDecoratedManifoldObjective}} (ρ_prime >= 0.25) && throw( ErrorException("ρ_prime must be strictly smaller than 0.25 but it is $ρ_prime.") ) @@ -474,8 +474,9 @@ function trust_regions!( augmentation_threshold=augmentation_threshold, (project!)=project!, ) - trs = decorate_state!(trs; kwargs...) - return get_solver_return(solve!(mp, trs)) + dtrs = decorate_state!(trs; kwargs...) + solve!(mp, dtrs) + return get_solver_return(get_objective(mp), dtrs) end function initialize_solver!(mp::AbstractManoptProblem, trs::TrustRegionsState) M = get_manifold(mp) diff --git a/test/functions/test_proximal_maps.jl b/test/functions/test_proximal_maps.jl index d540830a20..c44357437a 100644 --- a/test/functions/test_proximal_maps.jl +++ b/test/functions/test_proximal_maps.jl @@ -75,13 +75,13 @@ using Manifolds, Manopt, Test, Dates N2 = PowerManifold(M2, 3) pS, rS, qS = [-0.5, 0.1, 0.5] d = sum([pS, rS, qS] .* [1.0, -2.0, 1.0]) - m = min(0.3, abs(Manopt.sym_rem(d) / 6)) - s = sign(Manopt.sym_rem(d)) - pSc, rSc, qSc = Manopt.sym_rem.([pS, rS, qS] .- m .* s .* [1.0, -2.0, 1.0]) + m = min(0.3, abs(Manifolds.sym_rem(d) / 6)) + s = sign(Manifolds.sym_rem(d)) + pSc, rSc, qSc = Manifolds.sym_rem.([pS, rS, qS] .- m .* s .* [1.0, -2.0, 1.0]) pSr, rSr, qSr = prox_TV2(M2, 0.3, (pS, rS, qS)) @test sum(distance.(Ref(M2), [pSc, rSc, qSc], [pSr, rSr, qSr])) ≈ 0 # p=2 - t = 0.3 * Manopt.sym_rem(d) / (1 + 0.3 * 6.0) + t = 0.3 * Manifolds.sym_rem(d) / (1 + 0.3 * 6.0) @test sum( distance.( Ref(M2), diff --git a/test/plans/test_cache.jl b/test/plans/test_cache.jl new file mode 100644 index 0000000000..be859c0d6b --- /dev/null +++ b/test/plans/test_cache.jl @@ -0,0 +1,253 @@ +using LinearAlgebra, LRUCache, Manifolds, Manopt, Test + +include("../utils/dummy_types.jl") + +# Three dummy functors that are just meant to cound their calls +mutable struct TestCostCount + i::Int +end +TestCostCount() = TestCostCount(0) +function (tcc::TestCostCount)(M, p) + tcc.i += 1 + return norm(p) +end +mutable struct TestGradCount + i::Int +end +TestGradCount() = TestGradCount(0) +function (tgc::TestGradCount)(M, p) + tgc.i += 1 + return p +end +function (tgc::TestGradCount)(M, X, p) + tgc.i += 1 + X .= p + return X +end +mutable struct TestCostGradCount + i::Int +end +TestCostGradCount() = TestCostGradCount(0) +function (tcgc::TestCostGradCount)(M, p) + tcgc.i += 1 + return norm(p), p +end +function (tcgc::TestCostGradCount)(M, X, p) + tcgc.i += 1 + X .= p + return norm(p), X +end + +@testset "Test Caches" begin + @testset "Test Factory" begin + M = Euclidean(3) + # allocating + mgoa = ManifoldGradientObjective(TestCostCount(0), TestGradCount(0)) + s1 = objective_cache_factory(M, mgoa, :Simple) + @test s1 isa SimpleManifoldCachedObjective + @test objective_cache_factory(M, mgoa, :none) == mgoa + # pass a keyword + s2 = objective_cache_factory(M, mgoa, (:Simple, [], [:initialized => false])) + @test s2 isa SimpleManifoldCachedObjective + @test Manopt.is_objective_decorator(s2) + # but not initialized + @test !s2.X_valid + @test !s2.c_valid + # test fallbacks that do not decorate + @test objective_cache_factory(M, mgoa, :none) == mgoa + @test objective_cache_factory(M, mgoa, (:none, [])) == mgoa + @test objective_cache_factory(M, mgoa, (:none, [], [])) == mgoa + end + @testset "SimpleManifoldCachedObjective" begin + M = Euclidean(3) + p = zeros(3) + q = ones(3) + r = 2 * ones(3) + s = 3 * ones(3) + X = zero_vector(M, p) + # allocating + mgoa = ManifoldGradientObjective(TestCostCount(0), TestGradCount(0)) + mcgoa = ManifoldGradientObjective(TestCostCount(0), TestGradCount(0)) + sco1 = Manopt.SimpleManifoldCachedObjective(M, mgoa; p=p) + # We evaluated on init -> 1 + @test sco1.objective.gradient!!.i == 1 + @test sco1.objective.cost.i == 1 + @test get_gradient(M, sco1, p) == p + get_gradient!(M, X, sco1, p) + @test X == zero_vector(M, p) + @test get_cost(M, sco1, p) == norm(p) + # stil at 1 + @test sco1.objective.gradient!!.i == 1 + @test sco1.objective.cost.i == 1 + @test get_gradient(M, sco1, q) == q # triggers an evaluation + get_gradient!(M, X, sco1, q) # same point, copies + @test X == q + @test get_cost(M, sco1, q) == norm(q) + @test sco1.objective.cost.i == 2 + @test sco1.objective.gradient!!.i == 2 + # first grad! + get_gradient!(M, X, sco1, r) # triggers an evaluation + @test get_gradient(M, sco1, r) == X # cached + @test X == r + @test sco1.objective.gradient!!.i == 3 + @test get_cost_function(sco1) != get_cost_function(mgoa) + @test get_gradient_function(sco1) != get_gradient_function(mgoa) + + mgoi = ManifoldGradientObjective( + TestCostCount(0), TestGradCount(0); evaluation=InplaceEvaluation() + ) + sco2 = Manopt.SimpleManifoldCachedObjective(M, mgoi; p=p, initialized=false) + # We did not evaluate on init -> 1st eval + @test sco2.objective.gradient!!.i == 0 + @test sco2.objective.cost.i == 0 + @test get_gradient(M, sco2, p) == p + @test get_cost(M, sco2, p) == norm(p) + # now 1 + @test sco2.objective.gradient!!.i == 1 + @test sco2.objective.cost.i == 1 + # new point -> 2 + @test get_gradient(M, sco2, q) == q + get_gradient!(M, X, sco2, q) # cached + @test X == q + @test get_cost(M, sco2, q) == norm(q) + @test sco2.objective.gradient!!.i == 2 + @test sco2.objective.cost.i == 2 + get_gradient!(M, X, sco2, r) + @test get_gradient(M, sco2, r) == X # cached + @test X == r + + mcgoa = ManifoldCostGradientObjective(TestCostGradCount(0)) + sco3 = Manopt.SimpleManifoldCachedObjective(M, mcgoa; p=p, initialized=false) + # We do not evaluate on init -> still zero + @test sco3.objective.costgrad!!.i == 0 + @test get_gradient(M, sco3, p) == p + get_gradient!(M, X, sco3, p) + @test X == p + @test get_cost(M, sco3, p) == norm(p) + # stil at 1 + @test sco3.objective.costgrad!!.i == 1 + @test get_gradient(M, sco3, q) == q + get_gradient!(M, X, sco3, q) # cached + @test X == q + @test get_cost(M, sco3, q) == norm(q) # cached + @test sco3.objective.costgrad!!.i == 2 + get_gradient!(M, X, sco3, r) + @test X == r + @test get_gradient(M, sco3, r) == r # cached + @test get_cost(M, sco3, r) == norm(r) # cached + @test sco3.objective.costgrad!!.i == 3 + @test get_cost(M, sco3, s) == norm(s) + get_gradient!(M, X, sco3, s) # cached + @test X == s + @test get_gradient(M, sco3, s) == s # cached + @test sco3.objective.costgrad!!.i == 4 + + mcgoi = ManifoldCostGradientObjective( + TestCostGradCount(0); evaluation=InplaceEvaluation() + ) + sco4 = Manopt.SimpleManifoldCachedObjective(M, mcgoi; p=p) + # We evaluated on init -> evaluates twice + @test sco4.objective.costgrad!!.i == 2 + @test get_gradient(M, sco4, p) == p + get_gradient!(M, X, sco4, p) # cached + @test X == p + @test get_cost(M, sco4, p) == norm(p) + # stil at 2 + @test sco4.objective.costgrad!!.i == 2 + @test get_gradient(M, sco4, q) == q + get_gradient!(M, X, sco4, q) #cached + @test X == q + @test get_cost(M, sco4, q) == norm(q) + @test sco4.objective.costgrad!!.i == 3 + get_gradient!(M, X, sco4, r) + @test X == r + @test get_gradient(M, sco4, r) == r # cached + @test sco4.objective.costgrad!!.i == 4 + @test get_cost(M, sco4, s) == norm(s) + get_gradient!(M, X, sco4, s) # cached + @test X == s + @test get_gradient(M, sco4, s) == s # cached + @test sco4.objective.costgrad!!.i == 5 + end + @testset "ManifoldCachedObjective on Cost&Grad" begin + M = Sphere(2) + A = [2.0 1.0 0.0; 1.0 2.0 1.0; 0.0 1.0 2.0] + f(M, p) = p' * A * p + grad_f(M, p) = 2 * A * p + o = ManifoldGradientObjective(f, grad_f) + co = ManifoldCountObjective(M, o, [:Cost, :Gradient]) + lco = objective_cache_factory(M, co, (:LRU, [:Cost, :Gradient])) + ro = DummyDecoratedObjective(o) + #indecorated works as well + lco2 = objective_cache_factory(M, o, (:LRU, [:Cost, :Gradient])) + @test get_cost_function(lco2) != get_cost_function(o) + @test get_gradient_function(lco2) != get_gradient_function(o) + p = [1.0, 0.0, 0.0] + a = get_count(lco, :Cost) # usually 1 since creating lco calls that once + @test get_cost(M, lco, p) == 2.0 + @test get_cost(M, lco, p) == 2.0 + # but the second was cached so no cost eval + @test get_count(lco, :Cost) == a + 1 + # Gradient + b = get_count(lco, :Gradient) + X = get_gradient(M, lco, p) + @test X == grad_f(M, p) + # make sure this is safe, i.e. modifying X + X .= [1.0, 0.0, 1.0] + # does not affect the cache + @test get_gradient(M, lco, p) == grad_f(M, p) + X = get_gradient(M, lco, p) # restore X + Y = similar(X) + #Update Y inplace but without evaluating the gradient but taking it from the cache + get_gradient!(M, Y, lco, p) + @test Y == X + @test get_count(lco, :Gradient) == b + 1 + # + # CostGrad + f_f_grad(M, p) = (p' * A * p, 2 * A * p) + f_f_grad!(M, X, p) = (p' * A * p, X .= 2 * A * p) + o2a = ManifoldCostGradientObjective(f_f_grad) + co2a = ManifoldCountObjective(M, o2a, [:Cost, :Gradient]) + #pass size + lco2a = objective_cache_factory(M, co2a, (:LRU, [:Cost, :Gradient], 10)) + o2i = ManifoldCostGradientObjective(f_f_grad!; evaluation=InplaceEvaluation()) + co2i = ManifoldCountObjective(M, o2i, [:Cost, :Gradient]) + # pass keyword + lco2i = objective_cache_factory( + M, co2i, (:LRU, [:Cost, :Gradient], [:cache_size => 10]) + ) + # + c = get_count(lco2a, :Cost) # usually 1 since creating lco calls that once + @test get_cost(M, lco2a, p) == 2.0 + @test get_cost(M, lco2a, p) == 2.0 + # but the second was cached so no cost eval + @test get_count(lco2a, :Cost) == c + 1 + d = get_count(lco2a, :Gradient) + X = get_gradient(M, lco2a, p) + @test X == f_f_grad(M, p)[2] + Y = similar(X) + #Update Y inplace but without evaluating the gradient but taking it from the cache + get_gradient!(M, Y, lco, p) + @test Y == X + # But is Y also fixed in there ? not that we returned a reference to the cache. + Y .+= 1 + Z = similar(Y) + get_gradient!(M, Z, lco, p) + @test Z == X + get_gradient!(M, Y, lco, -p) #trigger cache with in-place + @test Y == -X + # Similar with + # Gradient cached already so no new evals + @test get_count(lco2a, :Gradient) == d + # Trigger caching on costgrad + X = get_gradient(M, lco2a, -p) + @test X == Y + # Trigger caching on costgrad! + get_gradient!(M, X, lco2i, -p) + @test X == Y + # Check default trigger + @test_throws DomainError Manopt.init_caches(M, [:Cost], Nothing) + @test_throws ErrorException Manopt.init_caches(M, [:None], LRU) + end + # Other tests are included with their respectives objective tests in the corresponding plans +end diff --git a/test/plans/test_cache_plan.jl b/test/plans/test_cache_plan.jl deleted file mode 100644 index bb829026c7..0000000000 --- a/test/plans/test_cache_plan.jl +++ /dev/null @@ -1,165 +0,0 @@ -using LinearAlgebra, Manifolds, Manopt, Test - -# Three dummy functors that are just meant to cound their calls -mutable struct TestCostCount - i::Int -end -TestCostCount() = TestCostCount(0) -function (tcc::TestCostCount)(M, p) - tcc.i += 1 - return norm(p) -end -mutable struct TestGradCount - i::Int -end -TestGradCount() = TestGradCount(0) -function (tgc::TestGradCount)(M, p) - tgc.i += 1 - return p -end -function (tgc::TestGradCount)(M, X, p) - tgc.i += 1 - X .= p - return X -end -mutable struct TestCostGradCount - i::Int -end -TestCostGradCount() = TestCostGradCount(0) -function (tcgc::TestCostGradCount)(M, p) - tcgc.i += 1 - return norm(p), p -end -function (tcgc::TestCostGradCount)(M, X, p) - tcgc.i += 1 - X .= p - return norm(p), X -end - -@testset "Test Caches" begin - @testset "Test Factory" begin - M = Euclidean(3) - # allocating - mgoa = ManifoldGradientObjective(TestCostCount(0), TestGradCount(0)) - s1 = objective_cache_factory(M, mgoa, :Simple) - @test s1 isa SimpleCacheObjective - @test objective_cache_factory(M, mgoa, :none) == mgoa - # pass a keyword - s2 = objective_cache_factory(M, mgoa, (:Simple, [], [:initialized => false])) - @test s2 isa SimpleCacheObjective - @test Manopt.is_objective_decorator(s2) - # but not initialized - @test !s2.X_valid - @test !s2.c_valid - @test objective_cache_factory(M, mgoa, (:none, [], [])) == mgoa - end - @testset "SimpleCacheObjective" begin - M = Euclidean(3) - p = zeros(3) - q = ones(3) - r = 2 * ones(3) - s = 3 * ones(3) - X = zero_vector(M, p) - # allocating - mgoa = ManifoldGradientObjective(TestCostCount(0), TestGradCount(0)) - mcgoa = ManifoldGradientObjective(TestCostCount(0), TestGradCount(0)) - sco1 = Manopt.SimpleCacheObjective(M, mgoa; p=p) - # We evaluated on init -> 1 - @test get_gradient_function(sco1).i == 1 - @test get_cost_function(sco1).i == 1 - @test get_gradient(M, sco1, p) == p - get_gradient!(M, X, sco1, p) - @test X == zero_vector(M, p) - @test get_cost(M, sco1, p) == norm(p) - # stil at 1 - @test get_gradient_function(sco1).i == 1 - @test get_cost_function(sco1).i == 1 - @test get_gradient(M, sco1, q) == q # triggers an evaluation - get_gradient!(M, X, sco1, q) # same point, copies - @test X == q - @test get_cost(M, sco1, q) == norm(q) - @test get_cost_function(sco1).i == 2 - @test get_gradient_function(sco1).i == 2 - # first grad! - get_gradient!(M, X, sco1, r) # triggers an evaluation - @test get_gradient(M, sco1, r) == X # cached - @test X == r - @test get_gradient_function(sco1).i == 3 - - mgoi = ManifoldGradientObjective( - TestCostCount(0), TestGradCount(0); evaluation=InplaceEvaluation() - ) - sco2 = Manopt.SimpleCacheObjective(M, mgoi; p=p, initialized=false) - # We did not evaluate on init -> 1st eval - @test get_gradient_function(sco2).i == 0 - @test get_cost_function(sco2).i == 0 - @test get_gradient(M, sco2, p) == p - @test get_cost(M, sco2, p) == norm(p) - # now 1 - @test get_gradient_function(sco2).i == 1 - @test get_cost_function(sco2).i == 1 - # new point -> 2 - @test get_gradient(M, sco2, q) == q - get_gradient!(M, X, sco2, q) # cached - @test X == q - @test get_cost(M, sco2, q) == norm(q) - @test get_gradient_function(sco2).i == 2 - @test get_cost_function(sco2).i == 2 - get_gradient!(M, X, sco2, r) - @test get_gradient(M, sco2, r) == X # cached - @test X == r - - mcgoa = ManifoldCostGradientObjective(TestCostGradCount(0)) - sco3 = Manopt.SimpleCacheObjective(M, mcgoa; p=p, initialized=false) - # We do not evaluate on init -> still zero - @test sco3.objective.costgrad!!.i == 0 - @test get_gradient(M, sco3, p) == p - get_gradient!(M, X, sco3, p) - @test X == p - @test get_cost(M, sco3, p) == norm(p) - # stil at 1 - @test sco3.objective.costgrad!!.i == 1 - @test get_gradient(M, sco3, q) == q - get_gradient!(M, X, sco3, q) # cached - @test X == q - @test get_cost(M, sco3, q) == norm(q) # cached - @test sco3.objective.costgrad!!.i == 2 - get_gradient!(M, X, sco3, r) - @test X == r - @test get_gradient(M, sco3, r) == r # cached - @test get_cost(M, sco3, r) == norm(r) # cached - @test sco3.objective.costgrad!!.i == 3 - @test get_cost(M, sco3, s) == norm(s) - get_gradient!(M, X, sco3, s) # cached - @test X == s - @test get_gradient(M, sco3, s) == s # cached - @test sco3.objective.costgrad!!.i == 4 - - mcgoi = ManifoldCostGradientObjective( - TestCostGradCount(0); evaluation=InplaceEvaluation() - ) - sco4 = Manopt.SimpleCacheObjective(M, mcgoi; p=p) - # We evaluated on init -> evaluates twice - @test sco4.objective.costgrad!!.i == 2 - @test get_gradient(M, sco4, p) == p - get_gradient!(M, X, sco4, p) # cached - @test X == p - @test get_cost(M, sco4, p) == norm(p) - # stil at 2 - @test sco4.objective.costgrad!!.i == 2 - @test get_gradient(M, sco4, q) == q - get_gradient!(M, X, sco4, q) #cached - @test X == q - @test get_cost(M, sco4, q) == norm(q) - @test sco4.objective.costgrad!!.i == 3 - get_gradient!(M, X, sco4, r) - @test X == r - @test get_gradient(M, sco4, r) == r # cached - @test sco4.objective.costgrad!!.i == 4 - @test get_cost(M, sco4, s) == norm(s) - get_gradient!(M, X, sco4, s) # cached - @test X == s - @test get_gradient(M, sco4, s) == s # cached - @test sco4.objective.costgrad!!.i == 5 - end -end diff --git a/test/plans/test_constrained_plan.jl b/test/plans/test_constrained_plan.jl index 29f850725a..a982a92a48 100644 --- a/test/plans/test_constrained_plan.jl +++ b/test/plans/test_constrained_plan.jl @@ -1,4 +1,6 @@ -using Manopt, ManifoldsBase, Test +using LRUCache, Manopt, ManifoldsBase, Test + +include("../utils/dummy_types.jl") @testset "Constrained Plan" begin M = ManifoldsBase.DefaultManifold(3) @@ -180,4 +182,207 @@ using Manopt, ManifoldsBase, Test end end end + @testset "Objetive Decorator passthrough" begin + for obj in [cofa, cofm, cova, covm] + ddo = DummyDecoratedObjective(obj) + @test get_constraints(M, ddo, p) == get_constraints(M, obj, p) + @test get_equality_constraints(M, ddo, p) == get_equality_constraints(M, obj, p) + @test get_inequality_constraints(M, ddo, p) == + get_inequality_constraints(M, obj, p) + Xe = get_grad_equality_constraints(M, ddo, p) + Ye = get_grad_equality_constraints(M, obj, p) + @test Ye == Xe + get_grad_equality_constraints!(M, Xe, ddo, p) + get_grad_equality_constraints!(M, Ye, obj, p) + @test Ye == Xe + for i in 1:1 #num of equ constr + @test get_equality_constraint(M, ddo, p, i) == + get_equality_constraint(M, obj, p, i) + X = get_grad_equality_constraint(M, ddo, p, i) + Y = get_grad_equality_constraint(M, obj, p, i) + @test X == Y + X = get_grad_equality_constraint!(M, X, ddo, p, i) + Y = get_grad_equality_constraint!(M, Y, obj, p, i) + @test X == Y + end + for j in 1:2 # num eq constr + @test get_inequality_constraint(M, ddo, p, j) == + get_inequality_constraint(M, obj, p, j) + X = get_grad_inequality_constraint(M, ddo, p, j) + Y = get_grad_inequality_constraint(M, obj, p, j) + @test X == Y + X = get_grad_inequality_constraint!(M, X, ddo, p, j) + Y = get_grad_inequality_constraint!(M, Y, obj, p, j) + @test X == Y + end + Xe = get_grad_inequality_constraints(M, ddo, p) + Ye = get_grad_inequality_constraints(M, obj, p) + @test Ye == Xe + get_grad_inequality_constraints!(M, Xe, ddo, p) + get_grad_inequality_constraints!(M, Ye, obj, p) + @test Ye == Xe + end + end + @testset "Count Objective" begin + ccofa = Manopt.objective_count_factory( + M, + cofa, + [ + :Constraints, + :InequalityConstraints, + :InequalityConstraint, + :EqualityConstraints, + :EqualityConstraint, + :GradInequalityConstraints, + :GradInequalityConstraint, + :GradEqualityConstraints, + :GradEqualityConstraint, + ], + ) + @test get_constraints(M, ccofa, p) == get_constraints(M, cofa, p) + @test get_count(ccofa, :Constraints) == 1 + @test get_equality_constraints(M, ccofa, p) == get_equality_constraints(M, cofa, p) + @test get_count(ccofa, :EqualityConstraints) == 1 + @test get_equality_constraint(M, ccofa, p, 1) == + get_equality_constraint(M, cofa, p, 1) + @test get_count(ccofa, :EqualityConstraint) == 1 + @test get_count(ccofa, :EqualityConstraint, 1) == 1 + @test get_inequality_constraints(M, ccofa, p) == + get_inequality_constraints(M, cofa, p) + @test get_count(ccofa, :InequalityConstraints) == 1 + @test get_inequality_constraint(M, ccofa, p, 1) == + get_inequality_constraint(M, cofa, p, 1) + @test get_inequality_constraint(M, ccofa, p, 2) == + get_inequality_constraint(M, cofa, p, 2) + @test get_count(ccofa, :InequalityConstraint) == [1, 1] + @test get_count(ccofa, :InequalityConstraint, 1) == 1 + @test get_count(ccofa, :InequalityConstraint, 2) == 1 + @test get_count(ccofa, :InequalityConstraint, [1, 2, 3]) == -1 + + Xe = get_grad_equality_constraints(M, cofa, p) + @test get_grad_equality_constraints(M, ccofa, p) == Xe + Ye = copy.(Ref(M), Ref(p), Xe) + get_grad_equality_constraints!(M, Ye, ccofa, p) + @test Ye == Xe + @test get_count(ccofa, :GradEqualityConstraints) == 2 + X = get_grad_equality_constraint(M, cofa, p, 1) + @test get_grad_equality_constraint(M, ccofa, p, 1) == X + Y = copy(M, p, X) + get_grad_equality_constraint!(M, Y, ccofa, p, 1) == X + @test Y == X + @test get_count(ccofa, :GradEqualityConstraint) == 2 + @test get_count(ccofa, :GradEqualityConstraint, 1) == 2 + Xi = get_grad_inequality_constraints(M, cofa, p) + @test get_grad_inequality_constraints(M, ccofa, p) == Xi + Yi = copy.(Ref(M), Ref(p), Xi) + @test get_grad_inequality_constraints!(M, Yi, ccofa, p) == Xi + @test get_count(ccofa, :GradInequalityConstraints) == 2 + X1 = get_grad_inequality_constraint(M, cofa, p, 1) + @test get_grad_inequality_constraint(M, ccofa, p, 1) == X1 + @test get_grad_inequality_constraint!(M, Y, ccofa, p, 1) == X1 + X2 = get_grad_inequality_constraint(M, cofa, p, 2) + @test get_grad_inequality_constraint(M, ccofa, p, 2) == X2 + @test get_grad_inequality_constraint!(M, Y, ccofa, p, 2) == X2 + @test get_count(ccofa, :GradInequalityConstraint) == [2, 2] + @test get_count(ccofa, :GradInequalityConstraint, 1) == 2 + @test get_count(ccofa, :GradInequalityConstraint, 2) == 2 + @test get_count(ccofa, :GradInequalityConstraint, [1, 2, 3]) == -1 + # test vectorial reset + reset_counters!(ccofa) + @test get_count(ccofa, :GradInequalityConstraint) == [0, 0] + end + @testset "Cache Objective" begin + cache_and_count = [ + :Constraints, + :InequalityConstraints, + :InequalityConstraint, + :EqualityConstraints, + :EqualityConstraint, + :GradInequalityConstraints, + :GradInequalityConstraint, + :GradEqualityConstraints, + :GradEqualityConstraint, + ] + ccofa = Manopt.objective_count_factory(M, cofa, cache_and_count) + cccofa = Manopt.objective_cache_factory(M, ccofa, (:LRU, cache_and_count)) + @test get_constraints(M, cofa, p) == get_constraints(M, cccofa, p) # counts + @test get_constraints(M, cofa, p) == get_constraints(M, cccofa, p) # cached + @test get_count(cccofa, :Constraints) == 1 + + ce = get_equality_constraints(M, cofa, p) + @test get_equality_constraints(M, cccofa, p) == ce # counts + @test get_equality_constraints(M, cccofa, p) == ce # cached + @test get_count(cccofa, :EqualityConstraints) == 1 + for i in 1 + ce_i = get_equality_constraint(M, cofa, p, i) + @test get_equality_constraint(M, cccofa, p, i) == ce_i # counts + @test get_equality_constraint(M, cccofa, p, i) == ce_i # cached + @test get_count(cccofa, :EqualityConstraint, i) == 1 + end + ci = get_inequality_constraints(M, cofa, p) + @test ci == get_inequality_constraints(M, cccofa, p) # counts + @test ci == get_inequality_constraints(M, cccofa, p) #cached + @test get_count(cccofa, :InequalityConstraints) == 1 + for j in 1:2 + ci_j = get_inequality_constraint(M, cofa, p, j) + @test get_inequality_constraint(M, cccofa, p, j) == ci_j # count + @test get_inequality_constraint(M, cccofa, p, j) == ci_j # cached + @test get_count(cccofa, :InequalityConstraint, j) == 1 + end + + Xe = get_grad_equality_constraints(M, cofa, p) + @test get_grad_equality_constraints(M, cccofa, p) == Xe # counts + @test get_grad_equality_constraints(M, cccofa, p) == Xe # cached + Ye = copy.(Ref(M), Ref(p), Xe) + get_grad_equality_constraints!(M, Ye, cccofa, p) # cached + @test Ye == Xe + @test get_count(ccofa, :GradEqualityConstraints) == 1 + Xe = get_grad_equality_constraints(M, cofa, -p) + get_grad_equality_constraints!(M, Ye, cccofa, -p) # ccounts + @test Ye == Xe + get_grad_equality_constraints!(M, Ye, cccofa, -p) # cached + @test Ye == Xe + @test get_grad_equality_constraints(M, cccofa, -p) == Xe # cached + for i in 1:1 + X = get_grad_equality_constraint(M, cofa, p, i) + @test get_grad_equality_constraint(M, cccofa, p, i) == X #counts + @test get_grad_equality_constraint(M, cccofa, p, i) == X #cached + Y = copy(M, p, X) + get_grad_equality_constraint!(M, Y, cccofa, p, i) == X # cached + @test Y == X + @test get_count(cccofa, :GradEqualityConstraint, i) == 1 + X = get_grad_equality_constraint(M, cofa, -p, i) + get_grad_equality_constraint!(M, Y, cccofa, -p, i) == X # counts + @test Y == X + get_grad_equality_constraint!(M, Y, cccofa, -p, i) == X # cached + @test Y == X + @test get_grad_equality_constraint(M, cccofa, -p, i) == X #cached + @test get_count(cccofa, :GradEqualityConstraint, i) == 2 + end + + Xi = get_grad_inequality_constraints(M, cofa, p) + @test get_grad_inequality_constraints(M, cccofa, p) == Xi # counts + @test get_grad_inequality_constraints(M, cccofa, p) == Xi # cached + Yi = copy.(Ref(M), Ref(p), Xi) + @test get_grad_inequality_constraints!(M, Yi, cccofa, p) == Xi # cached + @test get_count(cccofa, :GradInequalityConstraints) == 1 + Xi = get_grad_inequality_constraints(M, cofa, -p) + @test get_grad_inequality_constraints!(M, Yi, cccofa, -p) == Xi # counts + @test get_grad_inequality_constraints!(M, Yi, cccofa, -p) == Xi # cached + @test get_grad_inequality_constraints(M, cccofa, -p) == Xi # cached + @test get_count(cccofa, :GradInequalityConstraints) == 2 + for j in 1:2 + X = get_grad_inequality_constraint(M, cofa, p, j) + @test get_grad_inequality_constraint(M, cccofa, p, j) == X # counts + @test get_grad_inequality_constraint(M, cccofa, p, j) == X # cached + Y = copy(M, p, X) + @test get_grad_inequality_constraint!(M, Y, cccofa, p, j) == X # cached + @test get_count(ccofa, :GradInequalityConstraint, j) == 1 + X = get_grad_inequality_constraint(M, cofa, -p, j) + @test get_grad_inequality_constraint!(M, Y, cccofa, -p, j) == X # counts + @test get_grad_inequality_constraint!(M, Y, cccofa, -p, j) == X # cached + @test get_grad_inequality_constraint(M, cccofa, p, j) == X # cached + @test get_count(ccofa, :GradInequalityConstraint, j) == 2 + end + end end diff --git a/test/plans/test_counts.jl b/test/plans/test_counts.jl new file mode 100644 index 0000000000..b67e9c0c29 --- /dev/null +++ b/test/plans/test_counts.jl @@ -0,0 +1,56 @@ +using Manifolds, Manopt, Test + +include("../utils/dummy_types.jl") + +@testset "Counting Objective test" begin + M = Sphere(2) + A = [2.0 1.0 0.0; 1.0 2.0 1.0; 0.0 1.0 2.0] + f(M, p) = p' * A * p + grad_f(M, p) = project(M, p, 2 * A * p) + obj = ManifoldGradientObjective(f, grad_f) + c_obj = ManifoldCountObjective(M, obj, [:Cost, :Gradient]) + # function acessors are different since the right is still counting. + @test get_cost_function(obj) != get_cost_function(c_obj) + @test get_gradient_function(obj) != get_gradient_function(c_obj) + p = [1.0, 0.0, 0.0] + X = [1.0, 1.0, 0.0] + get_cost(M, c_obj, p) + @test get_count(c_obj, :Cost) == 1 + @test get_count(c_obj, :NonExistent) == -1 + Y = similar(X) + get_gradient(M, c_obj, p) + get_gradient!(M, Y, c_obj, p) + # both are counted + @test get_count(c_obj, :Gradient) == 2 + get_gradient(M, obj, p) + # others do not affect the counter + @test get_count(c_obj, :Gradient) == 2 + # also decorated objects can be wrapped to be counted + ro = DummyDecoratedObjective(obj) + c_obj2 = ManifoldCountObjective(M, ro, [:Gradient]) + get_gradient(M, c_obj2, p) + @test get_count(c_obj2, :Gradient) == 1 + @test_throws ErrorException get_count(ro, :Cost) # Errors since no CountObj + @test get_count(c_obj2, :Cost) == -1 # Does not count cost + @test_throws ErrorException get_count(c_obj2, :Cost, :error) + @test startswith(repr(c_obj), "## Statistics") + @test startswith(Manopt.status_summary(c_obj), "## Statistics") + # also if we get any (nonspecific tuple) - i.e. if the second is a point + @test startswith(repr((c_obj, p)), "## Statistics") + # but this also includes the hint, how to access the result + @test endswith(repr((c_obj, p)), "on this variable.") + rc_obj = DummyDecoratedObjective(c_obj) + @test get_count(rc_obj, :Gradient) == 2 #still works if count is encapsulated + @test_throws ErrorException get_count(obj, :Gradient) # no count objective + @test get_count(rc_obj, :Gradient, 1) == 2 #still works if count is encapsulated + @test_throws ErrorException get_count(obj, :Gradient, 1) # no count objective + # test fallbacks + @test get_count(c_obj, :None, 1) == -1 + @test get_count(c_obj, :Gradient, 2) == -1 # out of range + @test get_count(c_obj, :Gradient, [2, 1]) == -1 #nonfitting dimensions + reset_counters!(c_obj) + @test get_count(c_obj, :Gradient) == 0 + @test get_count(c_obj, :Cost) == 0 + reset_counters!(rc_obj) # also works on decorated counters + @test_throws ErrorException reset_counters!(obj) # errors on non-counter ones +end diff --git a/test/plans/test_difference_of_convex_plan.jl b/test/plans/test_difference_of_convex_plan.jl new file mode 100644 index 0000000000..189c0a28ec --- /dev/null +++ b/test/plans/test_difference_of_convex_plan.jl @@ -0,0 +1,62 @@ +using LRUCache, LinearAlgebra, Manopt, Manifolds, Test +include("../utils/dummy_types.jl") + +@testset "Difference of Convex Plan" begin + n = 2 + M = SymmetricPositiveDefinite(n) + g(M, p) = log(det(p))^4 + 1 / 4 + h(M, p) = log(det(p))^2 + grad_h(M, p) = 2 * log(det(p)) * p + f(M, p) = g(M, p) - h(M, p) + p = log(2) * Matrix{Float64}(I, n, n) + + dc_obj = ManifoldDifferenceOfConvexObjective(f, grad_h) + dcp_obj = ManifoldDifferenceOfConvexProximalObjective(grad_h; cost=f) + @testset "Objetive Decorator passthrough" begin + for obj in [dc_obj, dcp_obj] + ddo = DummyDecoratedObjective(obj) + X = get_subtrahend_gradient(M, ddo, p) + @test X == get_subtrahend_gradient(M, obj, p) + Y = zero_vector(M, p) + Z = zero_vector(M, p) + get_subtrahend_gradient!(M, Y, ddo, p) + get_subtrahend_gradient!(M, Z, obj, p) + @test Y == Z + end + end + @testset "Count" begin + for obj in [dc_obj, dcp_obj] + ddo = ManifoldCountObjective(M, obj, [:SubtrahendGradient]) + X = get_subtrahend_gradient(M, ddo, p) + @test X == get_subtrahend_gradient(M, obj, p) + Y = zero_vector(M, p) + Z = zero_vector(M, p) + get_subtrahend_gradient!(M, Y, ddo, p) + get_subtrahend_gradient!(M, Z, obj, p) + @test Y == Z + @test get_count(ddo, :SubtrahendGradient) == 2 + end + end + @testset "Cache" begin + for obj in [dc_obj, dcp_obj] + ddo = ManifoldCountObjective(M, obj, [:SubtrahendGradient]) + cddo = objective_cache_factory(M, ddo, (:LRU, [:SubtrahendGradient])) + X = get_subtrahend_gradient(M, obj, p) + @test X == get_subtrahend_gradient(M, cddo, p) + @test X == get_subtrahend_gradient(M, cddo, p) # Cached + Y = zero_vector(M, p) + get_subtrahend_gradient!(M, Y, cddo, p) # also cached + @test Y == X + @test get_count(ddo, :SubtrahendGradient) == 1 + + X = get_subtrahend_gradient(M, obj, 2 .* p) + Y = zero_vector(M, 2 .* p) + get_subtrahend_gradient!(M, Y, cddo, 2 .* p) + @test Y == X + get_subtrahend_gradient!(M, Y, cddo, 2 .* p) # cached + @test Y == X + @test X == get_subtrahend_gradient(M, cddo, 2 .* p) # Cached + @test get_count(ddo, :SubtrahendGradient) == 2 + end + end +end diff --git a/test/plans/test_gradient_plan.jl b/test/plans/test_gradient_plan.jl index c39cf92564..85dbd9d55b 100644 --- a/test/plans/test_gradient_plan.jl +++ b/test/plans/test_gradient_plan.jl @@ -1,4 +1,5 @@ using Manopt, ManifoldsBase, Test +include("../utils/dummy_types.jl") @testset "Gradient Plan" begin io = IOBuffer() @@ -13,7 +14,8 @@ using Manopt, ManifoldsBase, Test @test isapprox(M, p, get_gradient(gst), [1.0, 0.0]) f(M, q) = distance(M, q, p) .^ 2 grad_f(M, q) = -2 * log(M, q, p) - mp = DefaultManoptProblem(M, ManifoldGradientObjective(f, grad_f)) + mgo = ManifoldGradientObjective(f, grad_f) + mp = DefaultManoptProblem(M, mgo) @test get_initial_stepsize(mp, gst) == 1.0 @test get_stepsize(mp, gst, 1) == 1.0 @test get_last_stepsize(mp, gst, 1) == 1.0 @@ -73,5 +75,25 @@ using Manopt, ManifoldsBase, Test @test isapprox(M, p, X, get_gradient(M, mcgo!, p)) get_gradient!(M, Y, mcgo!, p) @test isapprox(M, p, X, Y) + cmcgo = ManifoldCountObjective(M, mcgo, [:Cost, :Gradient]) + @test get_cost(M, cmcgo, p) == get_cost(M, mcgo, p) + @test get_gradient(M, cmcgo, p) == get_gradient(M, mcgo, p) + get_gradient!(M, Y, cmcgo, p) + get_gradient!(M, X, mcgo, p) + # Now we called both 3 times + @test get_count(cmcgo, :Gradient) == 3 + @test get_count(cmcgo, :Cost) == 3 + end + @testset "Objetive Decorator passthrough" begin + ddo = DummyDecoratedObjective(mgo) + @test get_cost(M, mgo, p) == get_cost(M, ddo, p) + @test get_gradient(M, mgo, p) == get_gradient(M, ddo, p) + X = zero_vector(M, p) + Y = zero_vector(M, p) + get_gradient!(M, X, ddo, p) + get_gradient!(M, Y, ddo, p) + @test X == Y + @test get_gradient_function(ddo) == get_gradient_function(mgo) + @test get_cost_function(ddo) == get_cost_function(mgo) end end diff --git a/test/plans/test_hessian_plan.jl b/test/plans/test_hessian_plan.jl index b24fe1eca0..df2d470045 100644 --- a/test/plans/test_hessian_plan.jl +++ b/test/plans/test_hessian_plan.jl @@ -1,4 +1,5 @@ -using Manopt, Manifolds, Test +using LRUCache, Manopt, Manifolds, Test +include("../utils/dummy_types.jl") @testset "Hessian access functions" begin M = Euclidean(2) @@ -34,4 +35,72 @@ using Manopt, Manifolds, Test get_preconditioner!(mp, Y, p, X) @test Y == X end + @testset "Objetive Decorator passthrough" begin + Y1 = zero_vector(M, p) + Y2 = zero_vector(M, p) + for obj in [mho1, mho2, mho3, mho4] + ddo = DummyDecoratedObjective(obj) + @test get_hessian(M, obj, p, X) == get_hessian(M, ddo, p, X) + get_hessian!(M, Y1, obj, p, X) + get_hessian!(M, Y2, ddo, p, X) + @test Y1 == Y2 + @test get_preconditioner(M, obj, p, X) == get_preconditioner(M, ddo, p, X) + get_preconditioner!(M, Y1, obj, p, X) + get_preconditioner!(M, Y2, ddo, p, X) + @test Y1 == Y2 + end + end + @testset "Counting Objective" begin + Y1 = zero_vector(M, p) + Y2 = zero_vector(M, p) + for obj in [mho1, mho2, mho3, mho4] + cobj = Manopt.objective_count_factory(M, obj, [:Hessian, :Preconditioner]) + ddo = DummyDecoratedObjective(obj) + @test get_hessian(M, obj, p, X) == get_hessian(M, cobj, p, X) + get_hessian!(M, Y1, obj, p, X) + get_hessian!(M, Y2, cobj, p, X) + @test Y1 == Y2 + @test get_preconditioner(M, obj, p, X) == get_preconditioner(M, cobj, p, X) + get_preconditioner!(M, Y1, obj, p, X) + get_preconditioner!(M, Y2, cobj, p, X) + @test Y1 == Y2 + @test get_count(cobj, :Hessian) == 2 + @test get_count(cobj, :Preconditioner) == 2 + end + end + @testset "LRU Cache Objective" begin + Y = zero_vector(M, p) + for obj in [mho1, mho2, mho3, mho4] + cobj = Manopt.objective_count_factory(M, obj, [:Hessian, :Preconditioner]) + ccobj = Manopt.objective_cache_factory( + M, cobj, (:LRU, [:Hessian, :Preconditioner]) + ) + Z = get_hessian(M, obj, p, X) + @test get_hessian(M, ccobj, p, X) == Z + @test get_hessian(M, ccobj, p, X) == Z #cached + get_hessian!(M, Y, ccobj, p, X) #cached + @test Y == Z + @test get_count(ccobj, :Hessian) == 1 + Z = get_hessian(M, obj, -p, -X) + get_hessian!(M, Y, ccobj, -p, -X) #cached + @test Y == Z + @test get_hessian(M, ccobj, -p, -X) == Z #cached + @test get_count(ccobj, :Hessian) == 2 + + Z = get_preconditioner(M, obj, p, X) + @test get_preconditioner(M, ccobj, p, X) == Z + @test get_preconditioner(M, ccobj, p, X) == Z #cached + get_preconditioner!(M, Y, ccobj, p, X) #cached + @test Y == Z + @test get_count(ccobj, :Preconditioner) == 1 + Z = get_preconditioner(M, obj, -p, -X) + get_preconditioner!(M, Y, ccobj, -p, -X) + @test Y == Z + get_preconditioner!(M, Y, ccobj, -p, -X) # Cached + @test Y == Z + @test get_preconditioner(M, ccobj, -p, -X) == Z #cached + @test get_preconditioner(M, ccobj, -p, -X) == Z #cached + @test get_count(ccobj, :Preconditioner) == 2 + end + end end diff --git a/test/plans/test_objective.jl b/test/plans/test_objective.jl index 7967104434..0c33812ad9 100644 --- a/test/plans/test_objective.jl +++ b/test/plans/test_objective.jl @@ -1,16 +1,28 @@ using Manopt, Test -import Manopt: dispatch_objective_decorator -struct DummyDecoObjective{O<:AbstractManifoldObjective} <: - AbstractManifoldObjective{AllocatingEvaluation} - objective::O -end -dispatch_objective_decorator(::DummyDecoObjective) = Val(true) +include("../utils/dummy_types.jl") + +@testset "Objective" begin + @testset "Test decorator" begin + o = ManifoldCostObjective(x -> x) + d = DummyDecoratedObjective(o) + @test (get_objective(d) isa ManifoldCostObjective) + @test Manopt.is_objective_decorator(d) + @test !Manopt.is_objective_decorator(o) + end + @testset "ReturnObjective" begin + o = ManifoldCostObjective(x -> x) + r = Manopt.ReturnManifoldObjective(o) + @test repr(o) == "ManifoldCostObjective{AllocatingEvaluation}" + @test repr(r) == "ManifoldCostObjective{AllocatingEvaluation}" + @test Manopt.status_summary(o) == "ManifoldCostObjective{AllocatingEvaluation}" + @test Manopt.status_summary(r) == "ManifoldCostObjective{AllocatingEvaluation}" + @test repr((o, 1.0)) == """ + ManifoldCostObjective{AllocatingEvaluation} -@testset "Test decorator" begin - o = ManifoldCostObjective(x -> x) - d = DummyDecoObjective(o) - @test (get_objective(d) isa ManifoldCostObjective) - @test Manopt.is_objective_decorator(d) - @test !Manopt.is_objective_decorator(o) + To access the solver result, call `get_solver_result` on this variable.""" + d = DummyDecoratedObjective(o) + r2 = Manopt.ReturnManifoldObjective(d) + @test repr(r) == "ManifoldCostObjective{AllocatingEvaluation}" + end end diff --git a/test/plans/test_primal_dual_plan.jl b/test/plans/test_primal_dual_plan.jl index 959956e5be..c0c449e81b 100644 --- a/test/plans/test_primal_dual_plan.jl +++ b/test/plans/test_primal_dual_plan.jl @@ -1,5 +1,7 @@ using Manopt, Manifolds, ManifoldsBase, Test +include("../utils/dummy_types.jl") + @testset "Test primal dual plan" begin # # Perform an really easy test, just compute a mid point @@ -273,4 +275,41 @@ using Manopt, Manifolds, ManifoldsBase, Test @test length(get_record(r)) == 1 end end + @testset "Objetive Decorator passthrough" begin + # PD + pdmo = PrimalDualManifoldObjective( + f, prox_f, prox_g_dual, adjoint_DΛ; Λ=Λ, linearized_forward_operator=DΛ + ) + ro = DummyDecoratedObjective(pdmo) + q1 = get_primal_prox(M, ro, 0.1, p0) + q2 = get_primal_prox(M, pdmo, 0.1, p0) + @test q1 == q2 + get_primal_prox!(M, q1, ro, 0.1, p0) + get_primal_prox!(M, q2, pdmo, 0.1, p0) + @test q1 == q2 + Y1 = get_dual_prox(N, ro, n, 0.1, X0) + Y2 = get_dual_prox(N, pdmo, n, 0.1, X0) + @test Y1 == Y2 + get_dual_prox!(N, Y1, ro, n, 0.1, X0) + get_dual_prox!(N, Y2, pdmo, n, 0.1, X0) + @test Y1 == Y2 + Y1 = linearized_forward_operator(M, N, ro, m, p0, n) + Y2 = linearized_forward_operator(M, N, pdmol, m, p0, n) + @test Y1 == Y2 + linearized_forward_operator!(M, N, Y1, ro, m, p0, n) + linearized_forward_operator!(M, N, Y2, pdmol, m, p0, n) + @test Y1 == Y2 + Z1 = adjoint_linearized_operator(M, N, ro, m, n, X0) + Z2 = adjoint_linearized_operator(M, N, pdmol, m, n, X0) + @test Z1 == Z2 + adjoint_linearized_operator!(M, N, Z1, ro, m, n, X0) + adjoint_linearized_operator!(M, N, Z2, pdmol, m, n, X0) + @test Z1 == Z2 + s = forward_operator(M, N, ro, p0) + t = forward_operator(M, N, pdmo, p0) + @test s == t + forward_operator!(M, N, s, ro, p0) + forward_operator!(M, N, t, pdmo, p0) + @test Y1 == Y2 + end end diff --git a/test/plans/test_proximal_plan.jl b/test/plans/test_proximal_plan.jl new file mode 100644 index 0000000000..23313399a5 --- /dev/null +++ b/test/plans/test_proximal_plan.jl @@ -0,0 +1,67 @@ +using LRUCache, Manopt, Manifolds, Test +import Manopt: get_proximal_map, get_proximal_map! +function get_proximal_map(M, o::ManifoldProximalMapObjective, λ, p) + return get_proximal_map(M, o, λ, p, 1) +end +function get_proximal_map!(M, q, o::ManifoldProximalMapObjective, λ, p) + return get_proximal_map!(M, q, o, λ, p, 1) +end + +include("../utils/dummy_types.jl") + +@testset "Proximal Plan" begin + M = Euclidean(2) + p = [1.0, 2.0] + Q = [[2.0, 3.0], [3.0, 4.0]] + f(M, p) = sum(distance(M, p, q) for q in Q) + proxes_f = Tuple((N, λ, p) -> prox_distance(N, λ, q, p) for q in Q) + ppo = ManifoldProximalMapObjective(f, proxes_f) + @testset "Objetive Decorator passthrough" begin + dppo = DummyDecoratedObjective(ppo) + for i in 1:2 + @test get_proximal_map(M, ppo, 0.1, p, i) == + get_proximal_map(M, dppo, 0.1, p, i) + q = copy(M, p) + r = copy(M, p) + get_proximal_map!(M, q, ppo, 0.1, p, i) + get_proximal_map!(M, r, dppo, 0.1, p, i) + @test q == r + end + end + @testset "Counts" begin + cppo = ManifoldCountObjective(M, ppo, [:ProximalMap]) + q = get_proximal_map(M, cppo, 0.1, p, 1) + @test q == get_proximal_map(M, ppo, 0.1, p, 1) + q2 = copy(M, p) + get_proximal_map!(M, q2, cppo, 0.1, p, 1) + @test q2 == q + @test get_count(cppo, :ProximalMap, 1) == 2 + # the single ones have to be tricked a bit + cppo2 = ManifoldCountObjective(M, ppo, Dict([:ProximalMap => 0])) + @test q == get_proximal_map(M, cppo2, 0.1, p) + get_proximal_map!(M, q2, cppo2, 0.1, p) + @test q2 == q + @test get_count(cppo2, :ProximalMap) == 2 + end + @testset "Cache" begin + cppo = ManifoldCountObjective(M, ppo, [:ProximalMap]) + ccppo = objective_cache_factory(M, cppo, (:LRU, [:ProximalMap])) + for i in 1:2 + q = get_proximal_map(M, ppo, 0.1, p, i) + @test q == get_proximal_map(M, ccppo, 0.1, p, i) + @test q == get_proximal_map(M, ccppo, 0.1, p, i) # Cached + q2 = copy(M, p) + get_proximal_map!(M, q2, ccppo, 0.1, p, i) # Cached + @test q2 == q + @test get_count(ccppo, :ProximalMap, i) == 1 + + q = get_proximal_map(M, ppo, 0.2, -p, i) + get_proximal_map!(M, q2, ccppo, 0.2, -p, i) + @test q2 == q + get_proximal_map!(M, q2, ccppo, 0.2, -p, i) # Cached + @test q2 == q + @test q == get_proximal_map(M, ccppo, 0.2, -p, i) # Cached + @test get_count(ccppo, :ProximalMap, i) == 2 + end + end +end diff --git a/test/plans/test_state.jl b/test/plans/test_state.jl index e665218ddc..ee5a53660c 100644 --- a/test/plans/test_state.jl +++ b/test/plans/test_state.jl @@ -1,21 +1,17 @@ using Manifolds, Manopt, Test, ManifoldsBase -import Manopt: set_manopt_parameter! using Dates -struct TestStateProblem{M<:AbstractManifold} <: AbstractManoptProblem{M} end -mutable struct TestState <: AbstractManoptSolverState - storage::Vector{Float64} -end -TestState() = TestState([]) -set_manopt_parameter!(s::TestState, ::Val, v) = s +include("../utils/dummy_types.jl") + +struct NoIterateState <: AbstractManoptSolverState end @testset "Manopt Solver State" begin @testset "Generic State" begin M = Euclidean(3) - pr = TestStateProblem{typeof(M)}() - s = TestState() + pr = DummyStateProblem{typeof(M)}() + s = DummyState() @test repr(Manopt.ReturnSolverState(s)) == "ReturnSolverState($s)" - @test Manopt.status_summary(Manopt.ReturnSolverState(s)) == "TestState(Float64[])" + @test Manopt.status_summary(Manopt.ReturnSolverState(s)) == "DummyState(Float64[])" a = ArmijoLinesearch(M; initial_stepsize=1.0) @test get_last_stepsize(pr, s, a) == 1.0 @test get_initial_stepsize(a) == 1.0 @@ -28,13 +24,13 @@ set_manopt_parameter!(s::TestState, ::Val, v) = s ) @test get_initial_stepsize(dec_step) == 10.0 M = Euclidean(3) - pr = TestStateProblem{typeof(M)}() - @test dec_step(pr, TestState(), 1) == 10.0 - @test dec_step(pr, TestState(), 2) == 5.0 + pr = DummyStateProblem{typeof(M)}() + @test dec_step(pr, DummyState(), 1) == 10.0 + @test dec_step(pr, DummyState(), 2) == 5.0 end @testset "Decorator State" begin - s = TestState(zeros(3)) + s = DummyState(zeros(3)) r = RecordSolverState(s, RecordIteration()) d = DebugSolverState(s, DebugIteration()) ret = Manopt.ReturnSolverState(s) @@ -68,10 +64,12 @@ set_manopt_parameter!(s::TestState, ::Val, v) = s @test_throws ErrorException get_gradient(s) @test_throws ErrorException get_gradient(r) - @test_throws ErrorException get_iterate(s) - @test_throws ErrorException get_iterate(r) + @test isnan(get_iterate(s)) # dummy returns nan + @test isnan(get_iterate(r)) # dummy returns nan @test_throws ErrorException set_iterate!(s, M, 0) @test_throws ErrorException set_iterate!(r, M, 0) + s2 = NoIterateState() + @test_throws ErrorException get_iterate(s2) end @testset "Iteration and Gradient setters" begin M = Euclidean(3) @@ -90,4 +88,33 @@ set_manopt_parameter!(s::TestState, ::Val, v) = s set_gradient!(d2, M, p, X) @test d2.state.X == ones(3) end + @testset "Generic Objective and State solver returns" begin + f(M, p) = 1 + o = ManifoldCostObjective(f) + ro = Manopt.ReturnManifoldObjective(o) + ddo = DummyDecoratedObjective(o) + s = DummyState() + rs = Manopt.ReturnSolverState(s) + @test Manopt.get_solver_return(o, rs) == s #no ReturnObjective + # Return O & S + (a, b) = Manopt.get_solver_return(ro, rs) + @test a == o + @test b == s + # Return just S + Manopt.get_solver_return(ddo, rs) == s + # both as tuples and they return the iterate + @test isnan(get_solver_result((ro, rs))) + @test isnan(get_solver_result((o, rs))) + @test isnan(get_solver_result(ro, rs)) + @test isnan(get_solver_result(o, rs)) + # But also if the second is already some other type + @test isnan(get_solver_result((ro, NaN))) + @test isnan(get_solver_result((o, NaN))) + @test isnan(get_solver_result(ro, NaN)) + @test isnan(get_solver_result(o, NaN)) + # unless overwritten, objectives to not display in these tuples. + @test repr((o, s)) == repr(s) + # Passdown + @test repr((ro, s)) == repr(s) + end end diff --git a/test/plans/test_stochastic_gradient_plan.jl b/test/plans/test_stochastic_gradient_plan.jl new file mode 100644 index 0000000000..26be1fc732 --- /dev/null +++ b/test/plans/test_stochastic_gradient_plan.jl @@ -0,0 +1,100 @@ +using LinearAlgebra, LRUCache, Manopt, Manifolds, Test +include("../utils/dummy_types.jl") + +@testset "Stochastic Gradient Plan" begin + M = Sphere(2) + # 5 point mean + p = [0.0, 0.0, 1.0] + s = 1.0 + pts = [ + exp(M, p, X) for + X in [zeros(3), [s, 0.0, 0.0], [-s, 0.0, 0.0], [0.0, s, 0.0], [0.0, -s, 0.0]] + ] + f(M, y) = 1 / 2 * sum([distance(M, y, x)^2 for x in pts]) + sgrad_f1(M, y) = [-log(M, y, x) for x in pts] + sgrad_f2 = [((M, y) -> -log(M, y, x)) for x in pts] + msgo1 = ManifoldStochasticGradientObjective(sgrad_f1; cost=f) + msgo2 = ManifoldStochasticGradientObjective(sgrad_f2; cost=f) + @testset "Objetive Decorator passthrough" begin + X = zero_vector(M, p) + Y = zero_vector(M, p) + Xa = [zero_vector(M, p) for p in pts] + Ya = [zero_vector(M, p) for p in pts] + for obj in [msgo1, msgo2] + ddo = DummyDecoratedObjective(obj) + @test get_gradients(M, obj, p) == get_gradients(M, ddo, p) + get_gradients!(M, Xa, obj, p) + get_gradients!(M, Ya, ddo, p) + @test Xa == Ya + for i in 1:length(sgrad_f2) + @test get_gradient(M, obj, p, i) == get_gradient(M, ddo, p, i) + get_gradient!(M, X, obj, p, i) + get_gradient!(M, Y, ddo, p, i) + @test X == Y + end + end + end + @testset "Count Objetive" begin + X = zero_vector(M, p) + Y = zero_vector(M, p) + Xa = [zero_vector(M, p) for p in pts] + Ya = [zero_vector(M, p) for p in pts] + for obj in [msgo1, msgo2] + ddo = ManifoldCountObjective( + M, obj, [:StochasticGradient, :StochasticGradients] + ) + @test get_gradients(M, obj, p) == get_gradients(M, ddo, p) + get_gradients!(M, Xa, obj, p) + get_gradients!(M, Ya, ddo, p) + @test Xa == Ya + for i in 1:length(sgrad_f2) + @test get_gradient(M, obj, p, i) == get_gradient(M, ddo, p, i) + get_gradient!(M, X, obj, p, i) + get_gradient!(M, Y, ddo, p, i) + @test X == Y + @test get_count(ddo, :StochasticGradient, i) == 2 + end + @test get_count(ddo, :StochasticGradients) == 2 + end + end + @testset "Cache Objetive" begin + X = zero_vector(M, p) + Y = zero_vector(M, p) + Xa = [zero_vector(M, p) for p in pts] + Ya = [zero_vector(M, p) for p in pts] + for obj in [msgo1, msgo2] + ddo = ManifoldCountObjective( + M, obj, [:StochasticGradient, :StochasticGradients] + ) + cddo = objective_cache_factory( + M, ddo, (:LRU, [:StochasticGradient, :StochasticGradients]) + ) + Xa = get_gradients(M, obj, p) + @test Xa == get_gradients(M, cddo, p) # counts + @test Xa == get_gradients(M, cddo, p) # cached + get_gradients!(M, Ya, cddo, p) # cached + @test Xa == Ya + Xa = get_gradients(M, obj, -p) + get_gradients!(M, Ya, cddo, -p) # counts + @test Xa == Ya + get_gradients!(M, Ya, cddo, -p) # cached + @test Xa == Ya + @test Xa == get_gradients(M, cddo, -p) # cached + @test get_count(cddo, :StochasticGradients) == 2 + for i in 1:length(sgrad_f2) + X = get_gradient(M, obj, p, i) + @test X == get_gradient(M, cddo, p, i) # counts + @test X == get_gradient(M, cddo, p, i) # cached + get_gradient!(M, Y, cddo, p, i) # cached + @test X == Y + X = get_gradient(M, obj, -p, i) + get_gradient!(M, Y, cddo, -p, i) # counts + @test X == Y + get_gradient!(M, Y, cddo, -p, i) # cached + @test X == Y + @test X == get_gradient(M, cddo, -p, i) # cached + @test get_count(cddo, :StochasticGradient, i) == 2 + end + end + end +end diff --git a/test/plans/test_subgradient_plan.jl b/test/plans/test_subgradient_plan.jl new file mode 100644 index 0000000000..821c7da1c1 --- /dev/null +++ b/test/plans/test_subgradient_plan.jl @@ -0,0 +1,54 @@ +using LRUCache, Manopt, Manifolds, Test +include("../utils/dummy_types.jl") +@testset "Subgradient Plan" begin + M = Euclidean(2) + p = [1.0, 2.0] + f(M, q) = distance(M, q, p) + function ∂f(M, q) + if distance(M, p, q) == 0 + return zero_vector(M, q) + end + return -log(M, q, p) / max(10 * eps(Float64), distance(M, p, q)) + end + mso = ManifoldSubgradientObjective(f, ∂f) + @testset "Objetive Decorator passthrough" begin + ddo = DummyDecoratedObjective(mso) + @test get_cost(M, mso, p) == get_cost(M, ddo, p) + @test get_subgradient(M, mso, p) == get_subgradient(M, ddo, p) + X = zero_vector(M, p) + Y = zero_vector(M, p) + get_subgradient!(M, X, mso, p) + get_subgradient!(M, Y, ddo, p) + @test X == Y + end + @testset "Count" begin + ddo = ManifoldCountObjective(M, mso, [:SubGradient]) + @test get_subgradient(M, mso, p) == get_subgradient(M, ddo, p) + X = zero_vector(M, p) + Y = zero_vector(M, p) + get_subgradient!(M, X, mso, p) + get_subgradient!(M, Y, ddo, p) + @test X == Y + @test get_count(ddo, :SubGradient) == 2 + end + @testset "Cache" begin + ddo = ManifoldCountObjective(M, mso, [:SubGradient]) + cddo = objective_cache_factory(M, ddo, (:LRU, [:SubGradient])) + X = get_subgradient(M, mso, p) + @test get_subgradient(M, cddo, p) == X + @test get_subgradient(M, cddo, p) == X #Cached + Y = zero_vector(M, p) + get_subgradient!(M, Y, cddo, p) # Cached + @test X == Y + @test get_count(ddo, :SubGradient) == 1 + + X = get_subgradient(M, mso, -p) + Y = zero_vector(M, p) + get_subgradient!(M, Y, cddo, -p) + @test X == Y + get_subgradient!(M, Y, cddo, -p) # Cached + @test X == Y + @test get_subgradient(M, cddo, -p) == X #Cached + @test get_count(ddo, :SubGradient) == 2 + end +end diff --git a/test/runtests.jl b/test/runtests.jl index 4a253fa148..233e348ded 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,18 +8,23 @@ include("utils/example_tasks.jl") include("plans/test_problem.jl") include("plans/test_state.jl") include("plans/test_conjugate_gradient_plan.jl") + include("plans/test_counts.jl") include("plans/test_debug.jl") + include("plans/test_difference_of_convex_plan.jl") include("plans/test_storage.jl") - include("plans/test_cache_plan.jl") + include("plans/test_cache.jl") include("plans/test_nelder_mead_plan.jl") include("plans/test_gradient_plan.jl") include("plans/test_constrained_plan.jl") include("plans/test_hessian_plan.jl") include("plans/test_primal_dual_plan.jl") + include("plans/test_proximal_plan.jl") include("plans/test_higher_order_primal_dual_plan.jl") include("plans/test_record.jl") include("plans/test_stepsize.jl") + include("plans/test_stochastic_gradient_plan.jl") include("plans/test_stopping_criteria.jl") + include("plans/test_subgradient_plan.jl") end @testset "Function Tests " begin include("functions/test_adjoint_differentials.jl") diff --git a/test/solvers/test_gradient_descent.jl b/test/solvers/test_gradient_descent.jl index dc5f442c49..efb301ad35 100644 --- a/test/solvers/test_gradient_descent.jl +++ b/test/solvers/test_gradient_descent.jl @@ -145,6 +145,17 @@ using Manopt, Manifolds, Test r = gradient_descent!(M, f, grad_f, n5; return_state=true) @test isapprox(M, n5, n2) @test startswith(repr(r), "# Solver state for `Manopt.jl`s Gradient Descent") + # State and a count objective – putting stats behind print + n6 = gradient_descent( + M, + f, + grad_f, + pts[1]; + count=[:Gradient], + return_objective=true, + return_state=true, + ) + @test repr(n6) == "$(n6[2])\n\n$(n6[1])" end @testset "Warning when cost increases" begin M = Sphere(2) diff --git a/test/solvers/test_primal_dual_semismooth_Newton.jl b/test/solvers/test_primal_dual_semismooth_Newton.jl index 793f073963..af3835f538 100644 --- a/test/solvers/test_primal_dual_semismooth_Newton.jl +++ b/test/solvers/test_primal_dual_semismooth_Newton.jl @@ -188,4 +188,26 @@ using Manopt, Manifolds, ManifoldsBase, Test ) y2 = o2 @test x_hat ≈ y2 atol = 2 * 1e-7 + @testset "Objetive Decorator passthrough" begin + # PDNSSN additionals + pdmsno = PrimalDualManifoldSemismoothNewtonObjective( + f, prox_f, Dprox_F, prox_g_dual, Dprox_G_dual, DΛ, adjoint_DΛ; + ) + ro = DummyDecoratedObjective(pdmsno) + X = zero_vector(M, x0) + Y = get_differential_primal_prox(M, pdmsno, 0.1, x0, X) + Y2 = get_differential_primal_prox(M, ro, 0.1, x0, X) + @test Y == Y2 + get_differential_primal_prox!(M, Y, pdmsno, 0.1, x0, X) + get_differential_primal_prox!(M, Y2, ro, 0.1, x0, X) + @test Y == Y2 + + X = zero_vector(N, ξ0) + Y = get_differential_dual_prox(N, pdmsno, n, 0.1, ξ0, X) + Y2 = get_differential_dual_prox(N, ro, n, 0.1, ξ0, X) + @test Y == Y2 + get_differential_dual_prox!(N, Y, pdmsno, n, 0.1, ξ0, X) + get_differential_dual_prox!(N, Y2, ro, n, 0.1, ξ0, X) + @test Y == Y2 + end end diff --git a/test/utils/dummy_types.jl b/test/utils/dummy_types.jl new file mode 100644 index 0000000000..6ee0ec7e55 --- /dev/null +++ b/test/utils/dummy_types.jl @@ -0,0 +1,19 @@ +using Manopt +import Manopt: get_iterate, set_manopt_parameter! +struct DummyDecoratedObjective{E,O<:AbstractManifoldObjective} <: + Manopt.AbstractDecoratedManifoldObjective{E,O} + objective::O +end +function DummyDecoratedObjective( + o::O +) where {E<:AbstractEvaluationType,O<:AbstractManifoldObjective{E}} + return DummyDecoratedObjective{E,O}(o) +end + +struct DummyStateProblem{M<:AbstractManifold} <: AbstractManoptProblem{M} end +mutable struct DummyState <: AbstractManoptSolverState + storage::Vector{Float64} +end +DummyState() = DummyState([]) +get_iterate(::DummyState) = NaN +set_manopt_parameter!(s::DummyState, ::Val, v) = s diff --git a/tutorials/ConstrainedOptimization.qmd b/tutorials/ConstrainedOptimization.qmd index 97ef4616b5..7b2052ed6c 100644 --- a/tutorials/ConstrainedOptimization.qmd +++ b/tutorials/ConstrainedOptimization.qmd @@ -1,5 +1,6 @@ --- title: How to do Constrained Optimization +author: Ronny Bergmann --- This tutorial is a short introduction to using solvers for constraint optimisation in [`Manopt.jl`](https://manoptjl.org). diff --git a/tutorials/CountAndCache.qmd b/tutorials/CountAndCache.qmd new file mode 100644 index 0000000000..0864a66741 --- /dev/null +++ b/tutorials/CountAndCache.qmd @@ -0,0 +1,120 @@ +--- +title: "How to Count and Cache Function Calls" +author: Ronny Bergmann +--- + +In this tutorial, we want to investigate the caching and counting (i.e. statistics) features +of [Manopt.jl](https://manoptjl.org). We will reuse the optimization tasks from the +introductionary tutorial [Get Started: Optimize!](https://manoptjl.org/stable/tutorials/Optimize!.html). + +## Introduction + +There are surely many ways to keep track for example of how often the cost function is called, +for example with a [functor](https://docs.julialang.org/en/v1/manual/methods/#Function-like-objects), as we used in an example in [How to Record Data](https://manoptjl.org/stable/tutorials/HowtoRecord.html) + +```{julia} +#| eval: false +mutable struct MyCost{I<:Integer} + count::I +end +MyCost() = MyCost{Int64}(0) +function (c::MyCost)(M, x) + c.count += 1 + # [ .. Actual implementation of the cost here ] +end +``` + +This still leaves a bit of work to the user, especially for tracking more than just the number of cost function evaluations. + +When the a function like objective or gradient is expensive to compute, it may make sense to cache its results. +Manopt.jl tries to minimize the number of repeated calls but sometimes they are necessary and harmless when the function is cheap to compute. +Caching of expensive function calls can for example be added using [Memoize.jl](https://github.com/JuliaCollections/Memoize.jl) by the user. +The approach in the solvers of [Manopt.jl](https://manoptjl.org) aims to simplify adding +both these capabilities on the level of calling a solver. + +## Technical Background + +The two ingdredients for a solver in [Manopt.jl](https://manoptjl.org) +are the [`AbstractManoptProblem`](@ref) and the [`AbstractManoptSolverState`](@ref), where the +former consists of the domain, that is the [manifold](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/types.html#The-AbstractManifold) and [`AbstractManifoldObjective`](@ref). + +Both recording and debug capabilities are implemented in a decorator pattern to the solver state. +They can be easily added using the `record=` and `debug=` in any solver call. +This pattern was recently extended, such that also the objective can be decorated. +This is how both caching and counting are implemented, as decorators of the [`AbstractManifoldObjective`](@ref) +and hence for example changing/extending the behaviour of a call to [`get_cost`](@ref). + +Let's finish off the technical background by loading the necessary packages. +Besides [Manopt.jl](https://manoptjl.org) and [Manifolds.jl](https://juliamanifolds.github.io/Manifolds.jl/latest/) we also need +[LRUCaches.jl](https://github.com/JuliaCollections/LRUCache.jl) which are (since Julia 1.9) a weak dependency and provide +the _least recently used_ strategy for our caches. + +```{julia} +#| echo: false +#| code-fold: true +#| output: false +using Pkg; +cd(@__DIR__) +Pkg.activate("."); # for reproducibility use the local tutorial environment. +Pkg.develop(path="../") # a trick to work on the local dev version +``` + + +```{julia} +#| output: false +using Manopt, Manifolds, Random, LRUCache +``` + +## Counting + +We first define our task, the Riemannian Center of Mass from the [Get Started: Optimize!](https://manoptjl.org/stable/tutorials/Optimize!.html) tutorial. + +```{julia} +n = 100 +σ = π / 8 +M = Sphere(2) +p = 1 / sqrt(2) * [1.0, 0.0, 1.0] +data = [exp(M, p, σ * rand(M; vector_at=p)) for i in 1:n]; +f(M, p) = sum(1 / (2 * n) * distance.(Ref(M), Ref(p), data) .^ 2) +grad_f(M, p) = sum(1 / n * grad_distance.(Ref(M), data, Ref(p))); +``` + +to now count how often the cost and the gradient are called, we use the `count=` keyword +argument that works in any solver to specify the elements of the objective whose calls we +want to count calls to. A full list is available in the documentation of the +[`AbstractManifoldObjective`](@ref). +To also see the result, we have to set `return_objective=true`. +This returns `(objective, p)` instead of just the solver result `p`. +We can further also set `return_state=true` to get even more information about the solver run. + +```{julia} +gradient_descent(M, f, grad_f, data[1]; count=[:Cost, :Gradient], return_objective=true, return_state=true) +``` + +And we see that statistics are shown in the end. To now also cache these calls, +we can use the `cache=` keyword argument. +Since now both the cache and the count “extend” the functionality of the objective, +the order is important: On the high-level interface, the `count` is treated first, which +means that only actual function calls and not cache look-ups are counted. +With the proper initialisation, you can use any caches here that support the +`get!(function, cache, key)!` update. All parts of the objective that can currently be cached are listed at [`ManifoldCachedObjective`](@ref). The solver call has a keyword `cache` that takes a tuple`(c, vs, n)` of three arguments, where `c` is a symbol for the type of cache, `vs` is a vector of symbols, which calls to cache and `n` is the size of the cache. If the last element is not provided, a suitable default (currently`n=10`) is used. + +Here we want to use `c=:LRU` caches for `vs=[Cost, :Gradient]` with a size of `n=25`. + +```{julia} +r = gradient_descent(M, f, grad_f, data[1]; + count=[:Cost, :Gradient], + cache=(:LRU, [:Cost, :Gradient], 25), + return_objective=true, return_state=true) +``` + +Since the default setup with [`ArmijoLinesearch`](@ref) needs the gradient and the +cost, and similarly the stopping criterion might (independently) evaluate the gradient, +the caching is quite helpful here. + +And of course also for this advanced return value of the solver, we can still access the +result as usual: + +```{julia} +get_solver_result(r) +``` diff --git a/tutorials/GeodesicRegression.qmd b/tutorials/GeodesicRegression.qmd index 26a39fd421..76d5a1761d 100644 --- a/tutorials/GeodesicRegression.qmd +++ b/tutorials/GeodesicRegression.qmd @@ -1,5 +1,6 @@ --- title: How to perform Geodesic Regression +author: Ronny Bergmann --- Geodesic regression generalizes [linear regression](https://en.wikipedia.org/wiki/Linear_regression) diff --git a/tutorials/HowToRecord.qmd b/tutorials/HowToRecord.qmd index dd38f264df..5b176cf5f0 100644 --- a/tutorials/HowToRecord.qmd +++ b/tutorials/HowToRecord.qmd @@ -1,5 +1,6 @@ --- title: "How to Record Data During the Iterations" +author: Ronny Bergmann --- The recording and debugging features make it possible to record nearly any data during the iterations. diff --git a/tutorials/InplaceGradient.qmd b/tutorials/InplaceGradient.qmd index 25fe0e1a4b..9394c46c9b 100644 --- a/tutorials/InplaceGradient.qmd +++ b/tutorials/InplaceGradient.qmd @@ -1,5 +1,6 @@ --- title: "Speedup using Inplace Evaluation" +author: Ronny Bergmann --- When it comes to time critital operations, a main ingredient in Julia is given by diff --git a/tutorials/Optimize!.qmd b/tutorials/Optimize!.qmd index 9357efe1f1..7b86a63516 100644 --- a/tutorials/Optimize!.qmd +++ b/tutorials/Optimize!.qmd @@ -1,5 +1,6 @@ --- title: "Get Started: Optimize!" +author: Ronny Bergmann --- In this tutorial, we will both introduce the basics of optimisation on manifolds as well as diff --git a/tutorials/Project.toml b/tutorials/Project.toml index fb036d34af..6b90677548 100644 --- a/tutorials/Project.toml +++ b/tutorials/Project.toml @@ -2,6 +2,7 @@ BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000" +LRUCache = "8ac3fa9e-de4c-5943-b1dc-09c6b5f20637" ManifoldDiff = "af67fdf4-a580-4b9f-bbec-742ef357defd" Manifolds = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e" ManifoldsBase = "3362f125-f0bb-47a3-aa74-596ffd7ef2fb" @@ -10,8 +11,9 @@ Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" [compat] Colors = "0.12" -ManifoldDiff = "0.2" +ManifoldDiff = "0.3" Manifolds = "0.8.46" -ManifoldsBase = "0.13.30" +ManifoldsBase = "0.14.5" +LRUCache = "1.4" Manopt = "0.4.4" Plots = "1.38" diff --git a/tutorials/StochasticGradientDescent.qmd b/tutorials/StochasticGradientDescent.qmd index d0f83ee212..5e5bf302e3 100644 --- a/tutorials/StochasticGradientDescent.qmd +++ b/tutorials/StochasticGradientDescent.qmd @@ -1,5 +1,6 @@ --- title: How to Run Stochastic Gradient Descent +author: Ronny Bergmann --- This tutorial illustrates how to use the [`stochastic_gradient_descent`](https://manoptjl.org/stable/solvers/stochastic_gradient_descent.html) diff --git a/tutorials/_quarto.yml b/tutorials/_quarto.yml index ba1b778206..3c4294b2e5 100644 --- a/tutorials/_quarto.yml +++ b/tutorials/_quarto.yml @@ -22,4 +22,4 @@ format: variant: -raw_html+tex_math_dollars wrap: preserve -jupyter: julia-1.8 \ No newline at end of file +jupyter: julia-1.9 \ No newline at end of file