[WIP] Update parameters and implement grid search

Author: abelsiqueira

Project.toml

@@ -4,6 +4,7 @@ version = "0.1.0"
 
 [deps]
 ADNLPModels = "54578032-b7ea-4c30-94aa-7cbd1cce6c9a"
+Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 NLPModels = "a4795742-8479-5a88-8948-cc11e1c8c1a6"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 

solver-example.jl

@@ -21,4 +21,11 @@ output = solve!(solver, nlp, δ = 1e-2)
 
 #%%
 include("test/dummy_solver.jl")
-output, solver = DummySolver(nlp)
+output, solver = DummySolver(nlp)
+
+#%%
+include("test/dummy_solver.jl")
+problems = (
+  ADNLPModel(x -> (x[1] - a)^2 + b^2 * (x[2] - x[1]^2)^2, [-1.2; 1.0], x -> [x[1]^2 + x[2]^2 - 1], [0.0], [0.0]) for a = 0:0.1:1, b = 2:3
+)
+o = grid_search_tune(DummySolver, problems)

src/SolverCore.jl

@@ -1,12 +1,13 @@
 module SolverCore
 
 # stdlib
-using Printf
+using Logging, Printf
 
 # our packages
 using NLPModels
 
 include("solver.jl")
+include("grid-search-tuning.jl")
 include("logger.jl")
 include("stats.jl")
 

src/grid-search-tuning.jl

@@ -0,0 +1,99 @@
+export grid_search_tune
+
+"""
+    solver, results = grid_search_tune(SolverType, problems; kwargs...)
+
+Simple tuning of solver `SolverType` by grid search, on `problems`, which should be iterable.
+The following keyword arguments are available:
+- `success`: A function to be applied on a solver output that returns whether the problem has terminated succesfully. Defaults to `o -> o.status == :first_order`.
+- `costs`: A vector of cost functions and penalties. Each element is a tuple of two elements. The first is a function to be applied to the output of the solver, and the second is the cost when the solver fails (see `success` above) or throws an error. Defaults to
+```
+[
+  (o -> o.elapsed_time, 100.0),
+  (o -> o.counters.neval_obj + o.counters.neval_cons, 1000),
+  (o -> !success(o), 1),
+]
+```
+which represent the total elapsed_time (with a penalty of 100.0 for failures); the number of objective and constraints functions evaluations (with a penalty of 1000 for failures); and the number of failures.
+- `grid_length`: The number of points in the ranges of the grid for continuous points.
+- `solver_kwargs`: Arguments to be passed to the solver. Note: use this to set the stopping parameters, but not the other parameters being optimize.
+- Any parameters accepted by the `Solver`: a range to be used instead of the default range.
+
+The default ranges are based on the parameters types, and are as follows:
+- `:real`: linear range from `:min` to `:max` with `grid_length` points.
+- `:log`: logarithmic range from `:min` to `:max` with `grid_length` points. Computed by exp of linear range of `log(:min)` to `log(:max)`.
+- `:bool`: either `false` or `true`.
+- `:int`: integer range from `:min` to `:max`.
+"""
+function grid_search_tune(
+  ::Type{Solver},
+  problems;
+  success = o -> o.status == :first_order,
+  costs = [
+    (o -> o.elapsed_time, 100.0),
+    (o -> o.counters.neval_obj + o.counters.neval_cons, 1000),
+    (o -> !success(o), 1),
+  ],
+  grid_length = 10,
+  solver_kwargs = Dict(),
+  kwargs...
+) where Solver <: AbstractSolver
+
+  solver_params = parameters(Solver)
+  params = Dict()
+  for (k,v) in pairs(solver_params)
+    if v[:type] == :real
+      params[k] = LinRange(v[:min], v[:max], grid_length)
+    elseif v[:type] == :log
+      params[k] = exp.(LinRange(log(v[:min]), log(v[:max]), grid_length))
+    elseif v[:type] == :bool
+      params[k] = (false, true)
+    elseif v[:type] == :int
+      params[k] = v[:min]:v[:max]
+    end
+  end
+  for (k,v) in kwargs
+    params[k] = v
+  end
+
+  # Precompiling
+  nlp = first(problems)
+  try
+    solver = Solver(Val(:nosolve), nlp)
+    output = with_logger(NullLogger()) do
+      solve!(solver, nlp)
+    end
+  finally
+    finalize(nlp)
+  end
+
+  cost(θ) = begin
+    total_cost = [zero(x[2]) for x in costs]
+    for nlp in problems
+      reset!(nlp)
+      try
+        solver = Solver(Val(:nosolve), nlp)
+        output = with_logger(NullLogger()) do
+          solve!(solver, nlp; (k => θi for (k,θi) in zip(keys(solver_params), θ))...)
+        end
+        for (i, c) in enumerate(costs)
+          if success(output)
+            total_cost[i] += (c[1])(output)
+          else
+            total_cost[i] += c[2]
+          end
+        end
+      catch ex
+        for (i, c) in enumerate(costs)
+          total_cost[i] += c[2]
+        end
+        @show ex
+      finally
+        finalize(nlp)
+      end
+    end
+    total_cost
+  end
+
+  [θ => cost(θ) for θ in Iterators.product(values(params)...)]
+end

src/solver.jl

@@ -7,6 +7,10 @@ Base type for JSO-compliant solvers.
 """
 abstract type AbstractSolver{T} end
 
+function Base.show(io :: IO, solver :: AbstractSolver)
+    show(io, "Solver $(typeof(solver))")
+end
+
 """
     output = solve!(solver, problem)
 
@@ -27,7 +31,7 @@ Each key of `named_tuple` is the name of a parameter, and its value is a NamedTu
 - `default`: The default value of the parameter.
 - `type`: The type of the parameter, which can any of:
   - `:real`: A continuous value within a range
-  - `:log`: A continuous value that should be explorer logarithmically around it's lower value (usually 0) to avoid the bound itself.
+  - `:log`: A positive continuous value that should be explored logarithmically (like 10⁻², 10⁻¹, 1, 10).
   - `:int`: Integer value.
   - `:bool`: Boolean value.
 - `min`: Minimum value (may not be included for some parameter types).
@@ -36,32 +40,4 @@ Each key of `named_tuple` is the name of a parameter, and its value is a NamedTu
 function parameters(::Type{AbstractSolver{T}}) where T end
 
 parameters(::Type{S}) where S <: AbstractSolver = parameters(S{Float64})
-parameters(solver :: AbstractSolver) = parameters(typeof(solver))
-
-"""
-    nlp = parameter_problem(solver)
-
-Return the problem associated with the tuning of the parameters of `solver`.
-"""
-function parameter_problem(::AbstractSolver) end
-
-# parameter_problem(
-#   solver::DummySolver,
-#   problems,
-#   cost,
-#   cost_bad
-# ) = ADNLPModel(
-#   x -> begin
-#     total_cost = 0.0
-#     for nlp in problems
-#       try
-#         output = with_logger(NullLogger()) do
-#           output, _ = DummySolver(nlp)
-#         end
-#         total_cost += cost(output)
-#       catch
-#         total_cost +=
-#       end
-#     end
-#   end
-# )
+parameters(solver :: AbstractSolver) = parameters(typeof(solver))

test/dummy_solver.jl

@@ -1,5 +1,6 @@
 mutable struct DummySolver{T} <: AbstractSolver{T}
   initialized :: Bool
+  params :: Dict
   x :: Vector{T}
   xt :: Vector{T}
   gx :: Vector{T}
@@ -9,9 +10,31 @@ mutable struct DummySolver{T} <: AbstractSolver{T}
   ct :: Vector{T}
 end
 
-function DummySolver(::Type{T}, meta :: AbstractNLPModelMeta) where T
+function SolverCore.parameters(::Type{DummySolver{T}}) where T
+  (
+    α = (default=T(1e-2), type=:log, min=√√eps(T), max=one(T) / 2),
+    δ = (default=√eps(T), type=:log, min=√eps(T), max=√√√eps(T)),
+    reboot_y = (default=false, type=:bool)
+  )
+end
+
+function DummySolver(::Type{T}, meta :: AbstractNLPModelMeta; kwargs...) where T
   nvar, ncon = meta.nvar, meta.ncon
-  DummySolver{T}(true, zeros(T, nvar), zeros(T, nvar), zeros(T, nvar), zeros(T, nvar), zeros(T, ncon), zeros(T, ncon), zeros(T, ncon))
+  params = parameters(DummySolver{T})
+  solver = DummySolver{T}(true,
+    Dict(k => v[:default] for (k,v) in pairs(params)),
+    zeros(T, nvar),
+    zeros(T, nvar),
+    zeros(T, nvar),
+    zeros(T, nvar),
+    zeros(T, ncon),
+    zeros(T, ncon),
+    zeros(T, ncon),
+  )
+  for (k,v) in kwargs
+    solver.params[k] = v
+  end
+  solver
 end
 
 function DummySolver(::Type{T}, ::Val{:nosolve}, nlp :: AbstractNLPModel) where T
@@ -36,11 +59,16 @@ function SolverCore.solve!(solver::DummySolver{T}, nlp :: AbstractNLPModel;
   rtol :: Real = sqrt(eps(T)),
   max_eval :: Int = 1000,
   max_time :: Float64 = 30.0,
-  α :: Float64 = 1e-2,
-  δ :: Float64 = 1e-8,
+  kwargs...
 ) where T
   solver.initialized || error("Solver not initialized.")
   nvar, ncon = nlp.meta.nvar, nlp.meta.ncon
+  for (k,v) in kwargs
+    solver.params[k] = v
+  end
+  α = solver.params[:α]
+  δ = solver.params[:δ]
+  reboot_y = solver.params[:reboot_y]
 
   start_time = time()
   elapsed_time = 0.0
@@ -95,11 +123,16 @@ function SolverCore.solve!(solver::DummySolver{T}, nlp :: AbstractNLPModel;
     end
 
     x .= xt
-    y .+= t * Δy
+
 
     fx = ft
     grad!(nlp, x, gx)
     Jx = ncon > 0 ? jac(nlp, x) : zeros(T, 0, nvar)
+    if reboot_y
+      y .= -Jx' \ gx
+    else
+      y .+= t * Δy
+    end
     cx .= ct
     dual .= gx .+ Jx' * y
     elapsed_time = time() - start_time
@@ -132,10 +165,3 @@ function SolverCore.solve!(solver::DummySolver{T}, nlp :: AbstractNLPModel;
     iter=iter
   )
 end
-
-function SolverCore.parameters(::Type{DummySolver{T}}) where T
-  (
-    α = (default=T(1e-2), type=:log, min=zero(T), max=one(T)),
-    δ = (default=√eps(T), type=:log, min=zero(T), max=one(T)),
-  )
-end