Merge pull request #75 from JuliaNLSolvers/pkm/callbacks

Add callbacks throughout  and convergence test in neldermead

Author: pkofod

docs/mkdocs.yml

@@ -27,6 +27,7 @@ nav:
     - Tutorials:
         - Minimizing a function: 'optimization.md'
         - Solving non-linear equations: 'nonlineareq.md'
+        - Callbacks: 'callbacks.md'
 #        - Gradients and Hessians: 'user/gradientsandhessians.md'
 #        - Configurable Options: 'user/config.md'
 #        - Linesearch: 'algo/linesearch.md'

docs/src/callbacks.md

@@ -0,0 +1,161 @@
+# Callbacks
+
+NLSolvers.jl supports user-defined callbacks that are invoked after every iteration of any optimization solver. Callbacks let you monitor progress, log intermediate state, build a convergence trace, or stop the solver early based on custom criteria.
+
+## Basic usage
+
+A callback is any function (or callable object) that takes a single `info` argument and returns a `Bool`:
+
+```julia
+using NLSolvers
+
+callback = info -> begin
+    println("iter=$(info.iter), f=$(info.state.fz)")
+    return false   # continue optimization
+end
+
+solve(prob, x0, LineSearch(BFGS()),
+      OptimizationOptions(callback = callback))
+```
+
+Returning `true` stops the solver early; returning `false` lets it continue. The default `callback = nothing` disables the mechanism with zero overhead.
+
+## What the callback receives
+
+The `info` argument is a `NamedTuple` with three fields:
+
+| Field | Type | Description |
+|---|---|---|
+| `iter` | `Int` | Current iteration count (1-based) |
+| `time` | `Float64` | Elapsed seconds since `solve` started |
+| `state` | `NamedTuple` | Solver-specific state (see below) |
+
+## Solver-specific `state`
+
+The contents of `info.state` depend on which solver is running. Use `haskey` if you write generic callbacks.
+
+### Line search solvers (BFGS, DBFGS, DFP, SR1, L-BFGS, CG, Gradient Descent, Newton)
+
+`state` is the internal `objvars` named tuple:
+
+| Field | Description |
+|---|---|
+| `x` | Previous iterate |
+| `fx` | Objective value at `x` |
+| `∇fx` | Gradient at `x` |
+| `z` | New iterate (after line search) |
+| `fz` | Objective value at `z` |
+| `∇fz` | Gradient at `z` |
+| `B` | Current Hessian (or inverse) approximation; `nothing` for L-BFGS/CG |
+| `Pg` | Preconditioned gradient if using a preconditioner; otherwise `nothing` |
+
+### Trust region solvers
+
+Same fields as line search, plus:
+
+| Field | Description |
+|---|---|
+| `Δ` | Current trust region radius |
+| `rejected` | `true` if the previous step was rejected by the trust region rule |
+
+### Nelder-Mead
+
+| Field | Description |
+|---|---|
+| `simplex_vector` | Vector of simplex vertices |
+| `simplex_value` | Function values at each vertex |
+| `x_centroid` | Centroid of the simplex (excluding the worst vertex) |
+| `nm_obj` | Convergence metric — standard deviation of `simplex_value` |
+
+### Simulated Annealing
+
+| Field | Description |
+|---|---|
+| `x_best` | Best point found so far |
+| `f_best` | Best objective value found so far |
+| `x_now` | Current state of the chain |
+| `f_now` | Objective value at `x_now` |
+| `temperature` | Current temperature |
+
+### Particle Swarm
+
+| Field | Description |
+|---|---|
+| `X` | Current particle positions |
+| `X_best` | Each particle's personal best |
+| `Fs` | Function values at `X` |
+| `Fs_best` | Function values at `X_best` |
+| `x` | Global best particle |
+| `best_f` | Global best objective value |
+| `swarm_f` | Convergence metric for the swarm |
+
+### Brent's method (univariate)
+
+| Field | Description |
+|---|---|
+| `x` | Current best point |
+| `fx` | Function value at `x` |
+| `a`, `b` | Current bracketing interval `[a, b]` |
+| `v`, `w` | Two previous iterates |
+| `fv`, `fw` | Function values at `v` and `w` |
+
+### Active Box (projected Newton)
+
+| Field | Description |
+|---|---|
+| `x` | Previous iterate |
+| `z` | New iterate |
+| `fz` | Objective value at `z` |
+| `∇fz` | Gradient at `z` |
+| `B` | Hessian approximation |
+| `activeset` | Boolean vector indicating active bound constraints |
+
+## Examples
+
+### Build a convergence trace
+
+```julia
+trace = Float64[]
+solve(prob, x0, LineSearch(BFGS()),
+      OptimizationOptions(
+          callback = info -> (push!(trace, info.state.fz); false),
+          maxiter = 100,
+      ))
+```
+
+### Stop when the gradient is sufficiently small
+
+```julia
+gtol = 1e-6
+solve(prob, x0, LineSearch(BFGS()),
+      OptimizationOptions(
+          callback = info -> norm(info.state.∇fz, Inf) < gtol,
+          g_abstol = 0.0,  # disable built-in g-tolerance so callback wins
+      ))
+```
+
+### Time-limited optimization
+
+```julia
+time_limit = 5.0  # seconds
+solve(prob, x0, LineSearch(BFGS()),
+      OptimizationOptions(callback = info -> info.time > time_limit))
+```
+
+### Save iterates for plotting later
+
+Arrays in `info.state` (such as `z`, `∇fz`, `simplex_vector`) are aliases of live solver buffers — they will be overwritten on the next iteration. **Copy them if you need to retain them past the callback call.**
+
+```julia
+history = Vector{Vector{Float64}}()
+solve(prob, x0, LineSearch(BFGS()),
+      OptimizationOptions(
+          callback = info -> (push!(history, copy(info.state.z)); false),
+      ))
+```
+
+Scalar fields like `info.iter`, `info.time`, `info.state.fz` are values and do not need copying.
+
+## Performance
+
+The callback machinery has zero runtime overhead when `callback === nothing` (the default). The callback type is captured as a type parameter on `OptimizationOptions`, so the compiler eliminates the dispatch entirely. Passing a concrete callback function adds only the cost of calling that function once per iteration.

src/optimize/directsearch/neldermead.jl

@@ -193,7 +193,8 @@ function _solve(
     iter = 0
     nm_obj = f0
     is_converged = false
-    while iter <= options.maxiter && !any(is_converged)
+    callback_stopped = false
+    while iter <= options.maxiter && !any(is_converged) && !callback_stopped
         iter += 1
         nm_obj, x_centroid = iterate!(
             prob,
@@ -213,6 +214,7 @@ function _solve(
         if nm_obj ≤ options.nm_tol
             is_converged = true
         end
+        callback_stopped = _check_callback(options.callback, (iter=iter, time=time()-t0, state=(simplex_vector=simplex_vector, simplex_value=simplex_value, x_centroid=x_centroid, nm_obj=nm_obj)))
     end
     f_centroid_min = value(prob, x_centroid)
     f_min, i_f_min = findmin(simplex_value)
@@ -385,8 +387,9 @@ function _solve(
     is_converged = false
     iter = 0
     nm_obj = f0
+    callback_stopped = false
 
-    while iter <= options.maxiter && !any(is_converged)
+    while iter <= options.maxiter && !any(is_converged) && !callback_stopped
         iter += 1
 
         # Augment the iteration counter
@@ -467,10 +470,10 @@ function _solve(
         x_centroid = centroid(simplex_vector, i_order[end])
 
         nm_obj = nmobjective(simplex_value)
-        # if nm_x < 1e-18
-        #     break
-        # end
-        # check conv
+        if nm_obj ≤ options.nm_tol
+            is_converged = true
+        end
+        callback_stopped = _check_callback(options.callback, (iter=iter, time=time()-t0, state=(simplex_vector=simplex_vector, simplex_value=simplex_value, x_centroid=x_centroid, nm_obj=nm_obj)))
     end
     f_centroid_min = value(prob, x_centroid)
     f_min, i_f_min = findmin(simplex_value)

src/optimize/linesearch/conjugategradient.jl

@@ -236,13 +236,16 @@ function _solve(
     cgvars = CGVars(y, d, α, β, true)
 
     k = 1
+    callback_stopped = false
     objvars, P, cgvars = iterate(mstyle, cgvars, objvars, approach, problem, options)
     is_converged = converged(approach, objvars, ∇f0, options)
-    while k < options.maxiter && !any(is_converged)
+    callback_stopped = _check_callback(options.callback, (iter=k, time=time()-t0, state=objvars))
+    while k < options.maxiter && !any(is_converged) && !callback_stopped
         k += 1
         objvars, P, cgvars =
             iterate(mstyle, cgvars, objvars, approach, problem, options, P, false)
         is_converged = converged(approach, objvars, ∇f0, options)
+        callback_stopped = _check_callback(options.callback, (iter=k, time=time()-t0, state=objvars))
     end
     x, fx, ∇fx, z, fz, ∇fz, B = objvars
     return ConvergenceInfo(

src/optimize/linesearch/limitedquasinewton.jl

@@ -63,9 +63,11 @@ function _solve(
     objvars, qnvars, P =
         iterate(mstyle, 1, qnvars, objvars, P, approach, problem, problem, options)
     iter = 1
+    callback_stopped = false
     # Check for gradient convergence
     is_converged = converged(approach, objvars, ∇f0, options)
-    while iter <= options.maxiter && !any(is_converged)
+    callback_stopped = _check_callback(options.callback, (iter=iter, time=time()-t0, state=objvars))
+    while iter <= options.maxiter && !any(is_converged) && !callback_stopped
         iter += 1
         # take a step and update approximation
         objvars, qnvars, P = iterate(
@@ -82,6 +84,7 @@ function _solve(
         )
         # Check for gradient convergence
         is_converged = converged(approach, objvars, ∇f0, options)
+        callback_stopped = _check_callback(options.callback, (iter=iter, time=time()-t0, state=objvars))
     end
     x, fx, ∇fx, z, fz, ∇fz, B, Pg = objvars
     return ConvergenceInfo(

src/optimize/linesearch/quasinewton.jl

@@ -91,10 +91,12 @@ function _solve(
     ==============================#
     objvars, P, qnvars = iterate(mstyle, qnvars, objvars, P, approach, problem, options)
     iter = 1
+    callback_stopped = false
     # Check for gradient convergence
     is_converged = converged(approach, objvars, ∇f0, options)
     print_trace(approach, options, iter, t0, objvars)
-    while iter < options.maxiter && !any(is_converged)
+    callback_stopped = _check_callback(options.callback, (iter=iter, time=time()-t0, state=objvars))
+    while iter < options.maxiter && !any(is_converged) && !callback_stopped
         iter += 1
         #==============================
                      iterate
@@ -107,6 +109,7 @@ function _solve(
         ==============================#
         is_converged = converged(approach, objvars, ∇f0, options)
         print_trace(approach, options, iter, t0, objvars)
+        callback_stopped = _check_callback(options.callback, (iter=iter, time=time()-t0, state=objvars))
     end
     x, fx, ∇fx, z, fz, ∇fz, B, Pg = objvars
     return ConvergenceInfo(

src/optimize/problem_types.jl

@@ -199,7 +199,7 @@ function Base.show(io::IO, ci::ConvergenceInfo)
     println(io, "  Iterations:    $(info.iter)")
 end
 
-struct OptimizationOptions{T1,T2,T3,T4,Txn,Tgn}
+struct OptimizationOptions{T1,T2,T3,T4,Txn,Tgn,Tcb}
     x_abstol::T1
     x_reltol::T1
     x_norm::Txn
@@ -212,6 +212,7 @@ struct OptimizationOptions{T1,T2,T3,T4,Txn,Tgn}
     nm_tol::T3
     maxiter::T4
     show_trace::Bool
+    callback::Tcb
 end
 
 OptimizationOptions(;
@@ -222,11 +223,12 @@ OptimizationOptions(;
     g_reltol = 0.0,
     g_norm = x -> norm(x, Inf),
     f_limit = -Inf,
-    f_abstol = -Inf, # 
+    f_abstol = -Inf, #
     f_reltol = -Inf, # Not useful at 0 if for example we have quadric and trust region accept but objective is the same
     nm_tol = 1e-8,
     maxiter = 10000,
     show_trace = false,
+    callback = nothing,
 ) = OptimizationOptions(
     x_abstol,
     x_reltol,
@@ -240,8 +242,16 @@ OptimizationOptions(;
     nm_tol,
     maxiter,
     show_trace,
+    callback,
 )
 
+_check_callback(::Nothing, info) = false
+function _check_callback(cb, info)
+    result = cb(info)
+    result isa Bool || throw(ArgumentError("callback must return Bool, got $(typeof(result))"))
+    return result
+end
+
 struct MinResults{Tr,Tc<:ConvergenceInfo,Th,Ts,To}
     res::Tr
     conv::Tc

src/optimize/projected/projectednewton.jl

@@ -164,6 +164,9 @@ function _solve(
                 options,
             )
         end
+        if _check_callback(options.callback, (iter=iter, time=time()-t0, state=(x=x, z=z, fz=fz, ∇fz=∇fz, B=B, activeset=activeset)))
+            break
+        end
     end
     iter = options.maxiter
     return ConvergenceInfo(

src/optimize/randomsearch/particleswarm.jl

@@ -80,7 +80,8 @@ function solve(
     X_best[1] .= x0
     X[1] .= x0
     iter = 0
-    while iter < options.maxiter
+    callback_stopped = false
+    while iter < options.maxiter && !callback_stopped
         iter += 1
         limit_X!(X, lower, upper)
         Fs = batched_value(problem, Fs, X)
@@ -125,6 +126,7 @@ function solve(
         current_state, swarm_f = get_swarm_state(X, Fs, x, current_state)
         ω, c₁, c₂ = update_swarm_params!(c₁, c₂, ω, current_state, swarm_f)
         update_swarm!(X, X_best, x, n, V, ω, c₁, c₂)
+        callback_stopped = _check_callback(options.callback, (iter=iter, time=time()-t0, state=(X=X, X_best=X_best, Fs=Fs, Fs_best=Fs_best, x=x, best_f=best_f, swarm_f=swarm_f)))
     end
     best_f, x
     ConvergenceInfo(

src/optimize/randomsearch/simulatedannealing.jl

@@ -54,8 +54,9 @@ function solve(
     f0 = f_now
     temperature = f_best
     iter = 0
+    callback_stopped = false
     is_converged = converged(method, f_now, options)
-    while iter ≤ options.maxiter && !(is_converged)
+    while iter ≤ options.maxiter && !(is_converged) && !callback_stopped
         iter += 1
         # Determine the temperature for current iteration
         temperature = method.temperature(iter)
@@ -87,6 +88,7 @@ function solve(
             end
         end # n_per_temperature
         is_converged = converged(method, f_now, options)
+        callback_stopped = _check_callback(options.callback, (iter=iter, time=time()-t0, state=(x_best=x_best, f_best=f_best, x_now=x_now, f_now=f_now, temperature=temperature)))
     end
 
     ConvergenceInfo(