Allow prepending points. (#33)

Allow prepending points, fixes #31, #32.

Incidental: bump version

.github/workflows/CI.yml

@@ -17,7 +17,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.3' # Replace this with the minimum Julia version that your package supports. E.g. if your package requires Julia 1.5 or higher, change this to '1.5'.
+          - '1.6' # Replace this with the minimum Julia version that your package supports. E.g. if your package requires Julia 1.5 or higher, change this to '1.5'.
           - '1' # Leave this line unchanged. '1' will automatically expand to the latest stable 1.x release of Julia.
           - 'nightly'
         os:

Project.toml

@@ -1,7 +1,7 @@
 name = "MultistartOptimization"
 uuid = "3933049c-43be-478e-a8bb-6e0f7fd53575"
 authors = ["Tamas K. Papp <[email protected]>"]
-version = "0.1.3"
+version = "0.2.0"
 
 [deps]
 ArgCheck = "dce04be8-c92d-5529-be00-80e4d2c0e197"
@@ -13,12 +13,12 @@ Sobol = "ed01d8cd-4d21-5b2a-85b4-cc3bdc58bad4"
 
 [compat]
 ArgCheck = "1, 2"
-DocStringExtensions = "^0.8"
-NLopt = "0.5,0.6"
-Parameters = "^0.12"
+DocStringExtensions = "0.8"
+NLopt = "0.5, 0.6"
+Parameters = "0.12"
 Requires = "1.1"
-Sobol = "^1.3"
-julia = "^1.3"
+Sobol = "1.3"
+julia = "1.6"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

src/tiktak.jl

@@ -47,6 +47,29 @@ end
 """
 $(SIGNATURES)
 
+Determine if the value of the objective is a finite or `-Inf` real number.
+
+Internal, a basic sanity check to catch errors early.
+"""
+_acceptable_value(value) = value isa Real && (isfinite(value) || isinf(value)) # we don't want NaNs
+
+"""
+$(SIGNATURES)
+
+Evaluate `objective` at `location`, returning `location` and `value` in a `NamedTuple`.
+Perform basic sanity checks.
+
+Internal.
+"""
+function _objective_at_location(objective, location)
+        value = objective(location)
+        @argcheck _acceptable_value(value)
+        (; location, value)
+end
+
+"""
+$(SIGNATURES)
+
 Evaluate and return points of an `N`-element Sobol sequence.
 
 When `use_threads`, execution is parallelized using `Threads.@spawn`.
@@ -57,7 +80,7 @@ function sobol_starting_points(minimization_problem::MinimizationProblem, N::Int
     s = SobolSeq(lower_bounds, upper_bounds)
     skip(s, N)                  # better uniformity
     points = Iterators.take(s, N)
-    _initial(x) = (location = x, value = objective(x))
+    _initial(x) = _objective_at_location(objective, x)
     if use_threads
         map(fetch, map(x -> @spawn(_initial(x)), points))
     else
@@ -81,19 +104,28 @@ $(SIGNATURES)
 Solve `minimization_problem` by using `local_method` within `multistart_method`.
 
 When `use_threads`, initial point search is parallelized using `Threads.@spawn`.
+
+`prepend_points` should contain a vector of initial starting points that are prepended to
+the Sobol sequence. These are useful if a guess is available for the vicinity of the
+optimum.
 """
 function multistart_minimization(multistart_method::TikTak, local_method,
-                                 minimization_problem; use_threads = true)
+                                 minimization_problem;
+                                 use_threads = true,
+                                 prepend_points = Vector{Vector{Float64}}())
     @unpack quasirandom_N, initial_N, θ_min, θ_max, θ_pow = multistart_method
+    @unpack objective = minimization_problem
     quasirandom_points = sobol_starting_points(minimization_problem, quasirandom_N,
                                                use_threads)
     initial_points = _keep_lowest(quasirandom_points, initial_N)
+    all_points = vcat(map(x -> _objective_at_location(objective, x), prepend_points),
+                      initial_points)
     function _step(visited_minimum, (i, initial_point))
         θ = _weight_parameter(multistart_method, i)
         x = @. (1 - θ) * initial_point.location + θ * visited_minimum.location
         local_minimum = local_minimization(local_method, minimization_problem, x)
         local_minimum ≡ nothing && return visited_minimum
         local_minimum.value < visited_minimum.value ? local_minimum : visited_minimum
     end
-    foldl(_step, enumerate(initial_points); init = first(initial_points))
+    foldl(_step, enumerate(Iterators.drop(all_points, 1)); init = first(all_points))
 end

test/runtests.jl

@@ -43,15 +43,15 @@ end
 @testset "local gradient-based methods" begin
     function autodiff(fn)
         # adapted from here https://github.com/JuliaOpt/NLopt.jl/issues/128
-        function f(x) 
-            return fn(x) 
+        function f(x)
+            return fn(x)
         end
-    
+
         function f(x,∇f)
             if !(∇f == nothing) && (length(∇f) != 0)
                 ForwardDiff.gradient!(∇f,fn,x)
             end
-        
+
             fx = fn(x)
             return fx
         end
@@ -70,3 +70,23 @@ end
     end
 
 end
+
+@testset "prepending points and sanity checks" begin
+    N = 3
+    lb = .-ones(N)
+    ub = 2 .* ones(N)
+    vz1 = -5.0
+    z1 = lb .+ 1/√2 .* (ub .- lb) # special-cased, will be vz1
+    g = function(x)
+        # very narrow optimum near z1
+        α = 1 - exp(-sum(abs2, 100 .* (x .- z1)))
+        α * sum(abs2, x) + (1 - α) * vz1
+    end
+    P = MinimizationProblem(g, lb, ub)
+    MM, LM = TikTak(100), NLoptLocalMethod(NLopt.LN_BOBYQA)
+    r0 = multistart_minimization(MM, LM, P; use_threads = false)
+    @test r0.value ≈ 0 atol = 1e-9         # sanity check
+    r1 = multistart_minimization(MM, LM, P; use_threads = false, prepend_points = [z1])
+    @test r1.value == vz1
+    @test r1.location == z1
+end