Rel 3.3.0

* Add maximum time termination to options

Project.toml

@@ -1,7 +1,7 @@
 name = "NLLSsolver"
 uuid = "50b5cf9e-bf7a-4e29-8981-4280cbba70cb"
 authors = ["Oliver Woodford"]
-version = "3.2.3"
+version = "3.3.0"
 
 [deps]
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"

src/linearsolver.jl

@@ -30,3 +30,7 @@ solve!(ls::MultiVariateLSsparse, options) = ldiv!(ls.x, ldl_factorize!(ls.hessia
 solve!(ls::MultiVariateLSdense, options) = solve!(ls.x, ls.A.data, ls.b)
 solve!(x, A, b) = try_cholesky!(x, A, b)
 solve!(x, A::SparseMatrixCSC, b) = ldiv!(x, ldl(A), b)
+
+# Symmetric matrix inversion
+invsym(A::AbstractMatrix) = inv(bunchkaufman(A))
+invsym(A::StaticMatrix) = Size(A)[1] <= 14 ? inv(A) : inv(bunchkaufman(A))

src/marginalize.jl

@@ -10,7 +10,7 @@ function marginalize!(to::MultiVariateLS, from::MultiVariateLSsparse, blockind::
     N = length(blocks) - 1
     dataindices = view(from.A.indicestransposed.nzval, ind)
     # Compute inverse for the diagonal block (to be marginalized)
-    inverseblock = inv(reshape(view(from.A.data, (0:blocksz*blocksz-1) .+ dataindices[end]), blocksz, blocksz))
+    inverseblock = invsym(reshape(view(from.A.data, (0:blocksz*blocksz-1) .+ dataindices[end]), blocksz, blocksz))
     # For each non-marginalized block
     blockgrad = view(from.b, (0:blocksz-1) .+ from.boffsets[blockind])
     for a in 1:N
@@ -40,7 +40,7 @@ function marginalize!(to::MultiVariateLS, from::MultiVariateLSsparse, blockind::
     N = length(blocks) - 1
     dataindices = view(from.A.indicestransposed.nzval, ind)
     # Compute inverse for the diagonal block (to be marginalized)
-    inverseblock = inv(SizedMatrix{blocksz, blocksz}(view(from.A.data, SR(0, blocksz*blocksz-1).+dataindices[end])))
+    inverseblock = invsym(SizedMatrix{blocksz, blocksz}(view(from.A.data, SR(0, blocksz*blocksz-1).+dataindices[end])))
     # For each non-marginalized block
     blockgrad = SizedVector{blocksz}(view(from.b, SR(0, blocksz-1) .+ from.boffsets[blockind]))
     for a in 1:N

src/optimize.jl

@@ -62,6 +62,7 @@ function optimizeinternal!(problem::NLLSProblem, options::NLLSOptions, data, ite
     if length(problem.variables) != length(problem.varnext)
         problem.varnext = copy(problem.variables)
     end
+    stoptime = starttime + options.maxtime
     data.timeinit += Base.time_ns() - starttime
     # Initialize the linear problem
     data.timegradient += @elapsed_ns data.bestcost = costgradhess!(data.linsystem, problem.variables, problem.costs)
@@ -117,7 +118,8 @@ function optimizeinternal!(problem::NLLSProblem, options::NLLSOptions, data, ite
         converged |= (maxstep < options.dstep)                       << 6 # Max of the step size is too small
         converged |= (fails > options.maxfails)                      << 7 # Max number of consecutive failed iterations reach
         converged |= (data.iternum >= options.maxiters)              << 8 # Max number of iterations reached
-        converged |= terminate                                       << 9 # Terminated by the user-defined callback
+        converged |= (Base.time_ns() > stoptime)                     << 9 # Max amount of time exceeded
+        converged |= terminate                                       << 16 # Terminated by the user-defined callback (room left for new flags above)
         data.converged = converged
         if converged != 0
             break

src/structs.jl

@@ -25,12 +25,13 @@ struct NLLSOptions{T}
     dstep::Float64              # Minimum L-infinity norm of the update vector required to avoid termination
     maxfails::Int               # Maximum number of consecutive iterations that have a higher cost than the current best before termination
     maxiters::Int               # Maximum number of outer iterations
+    maxtime::UInt64             # Maximum optimization time allowed, in nano-seconds (converted from seconds in the constructor)
     iterator::NLLSIterator      # Inner iterator (see above for options)
     callback::T                 # Callback called every outer iteration - (cost, problem, data) -> (newcost, terminate::Bool) where terminate == true ends the optimization
     storecosts::Bool            # Indicates whether the cost per outer iteration should be stored
     storetrajectory::Bool       # Indicates whether the step per outer iteration should be stored
 end
-function NLLSOptions(; maxiters=100, reldcost=1.e-15, absdcost=1.e-15, dstep=1.e-15, maxfails=3, iterator=levenbergmarquardt, callback::T=nothing, storecosts=false, storetrajectory=false) where T
+function NLLSOptions(; maxiters=100, reldcost=1.e-15, absdcost=1.e-15, dstep=1.e-15, maxfails=3, maxtime=30.0, iterator=levenbergmarquardt, callback::T=nothing, storecosts=false, storetrajectory=false) where T
     if iterator == gaussnewton
         Base.depwarn("gaussnewton is deprecated. Use newton instead", :NLLSOptions)
     end
@@ -40,7 +41,7 @@ function NLLSOptions(; maxiters=100, reldcost=1.e-15, absdcost=1.e-15, dstep=1.e
     if !isnothing(callback)
         Base.depwarn("Setting callback in options is deprecated. Pass the callback directly to optimize!() instead", :NLLSOptions)
     end
-    NLLSOptions{T}(reldcost, absdcost, dstep, maxfails, maxiters, iterator, callback, storecosts, storetrajectory)
+    NLLSOptions{T}(reldcost, absdcost, dstep, maxfails, maxiters, UInt64(round(maxtime * 1e9)), iterator, callback, storecosts, storetrajectory)
 end
 
 struct NLLSResult
@@ -51,7 +52,7 @@ struct NLLSResult
     timecost::Float64                       # Time (in seconds) spent computing the cost
     timegradient::Float64                   # Time (in seconds) spent computing the residual gradients and constructing the linear problems
     timesolver::Float64                     # Time (in seconds) spent solving the linear problems
-    termination::Int                        # Set of flags indicating which termination criteria were met
+    termination::Int                        # Set of flags indicating which termination criteria were met - the value should not be relied upon
     niterations::Int                        # Number of outer optimization iterations performed
     costcomputations::Int                   # Number of cost computations performed
     gradientcomputations::Int               # Number of residual gradient computations performed
@@ -84,7 +85,8 @@ function Base.show(io::IO, x::NLLSResult)
     if 0 != x.termination & (1 << 6); println(io, "   Step size below threshold."); end
     if 0 != x.termination & (1 << 7); println(io, "   Too many consecutive iterations increasing the cost."); end
     if 0 != x.termination & (1 << 8); println(io, "   Maximum number of outer iterations reached."); end
-    userflags = x.termination >> 9
+    if 0 != x.termination & (1 << 9); println(io, "   Maximum allowed computation time exceeded."); end
+    userflags = x.termination >> 16
     if 0 != userflags; println(io, "   Terminated by user-defined callback, with flags: ", string(userflags, base=2)); end
 end
 

test/functional.jl

@@ -47,9 +47,9 @@ Base.eltype(::RosenbrockB) = Float64
     @test NLLSsolver.countcosts(NLLSsolver.resnum, subprob.costs) == 1
     @test NLLSsolver.cost(subprob) == 0.
 
-    # Check callback termination
-    result = NLLSsolver.optimize!(problem, NLLSsolver.NLLSOptions(), nothing, (cost, unusedargs...)->(cost, 13))
-    @test result.termination == (13 << 9)
+    # Check callback and max time termination
+    result = NLLSsolver.optimize!(problem, NLLSsolver.NLLSOptions(maxtime=0.0), nothing, (cost, unusedargs...)->(cost, 13))
+    @test result.termination == (1 << 9) | (13 << 16)
     @test result.niterations == 1
 
     # Optimize using Newton