Heuristics Interface (#162)

examples/nonlinear.jl

@@ -3,6 +3,7 @@ using LinearAlgebra
 import MathOptInterface
 using Random
 using Boscia
+using Bonobo
 import Bonobo
 using Printf
 using Dates
@@ -66,6 +67,52 @@ end
 # benchmarking Oracles
 FrankWolfe.benchmark_oracles(f, grad!, () -> rand(n), lmo; k=100)
 
-x, _, _ = Boscia.solve(f, grad!, lmo, verbose=true, print_iter=500)
+#################
+
+
+# are these lmo calls counted as well?
+
+# #####
+# # follow the gradient for a fixed number of steps and collect solutions on the way
+# #####
+
+# function follow_gradient_heuristic(tree::Bonobo.BnBTree, blmo::Boscia.BoundedLinearMinimizationOracle, x, k)
+#     nabla = similar(x)
+#     x_new = copy(x)
+#     sols = []
+#     for i in 1:k
+#         tree.root.problem.g(nabla,x_new)
+#         x_new = Boscia.compute_extreme_point(blmo, nabla)
+#         push!(sols, x_new)
+#     end
+#     return sols, false
+# end
+
+# #####
+# # rounding respecting the hidden feasible region structure
+# #####
+
+# function rounding_lmo_01_heuristic(tree::Bonobo.BnBTree, blmo::Boscia.BoundedLinearMinimizationOracle, x)
+#     nabla = zeros(length(x))
+#     for idx in tree.branching_indices
+#         nabla[idx] = 1 - 2*round(x[idx]) # (0.7, 0.3) -> (1, 0) -> (-1, 1) -> min -> (1,0)
+#     end
+#     x_rounded = Boscia.compute_extreme_point(blmo, nabla)
+#     return [x_rounded], false
+# end
+
+#####
+# geometric scaling like for a couple of steps
+#####
+
+
+depth = 5
+heu  = Boscia.Heuristic((tree, blmo, x) -> Boscia.follow_gradient_heuristic(tree,blmo,x, depth), 0.8, :follow_gradient)
+heu2 = Boscia.Heuristic(Boscia.rounding_lmo_01_heuristic, 0.8, :lmo_rounding)
+
+heuristics = [heu, heu2]
+# heuristics = []
+
+x, _, _ = Boscia.solve(f, grad!, lmo, verbose=true, print_iter=500, custom_heuristics=heuristics)
 
 @show x

examples/portfolio.jl

@@ -4,10 +4,12 @@ using Test
 using Random
 using SCIP
 using LinearAlgebra
-using Distributions
 import MathOptInterface
 const MOI = MathOptInterface
 
+# seed = 0x946d4b7835e92ffa takes 90 minutes to solve! -> not anymore
+seed = 0x946d4b7835e92ffa
+Random.seed!(seed)
 
 n = 30
 const ri = rand(n)
@@ -53,10 +55,12 @@ const Mi = (Ai + Ai') / 2
         return storage
     end
 
-    x, _, result = Boscia.solve(f, grad!, lmo, verbose=true, time_limit=600)
+    depth = 5
+    heu  = Boscia.Heuristic((tree, blmo, x) -> Boscia.follow_gradient_heuristic(tree,blmo,x, depth), 0.2, :follow_gradient)
+    heuristics = [heu]
+    # heuristics = []
+
+    x, _, result = Boscia.solve(f, grad!, lmo, verbose=true, time_limit=600, custom_heuristics=heuristics)
     @test dot(ai, x) <= bi + 1e-2
     @test f(x) <= f(result[:raw_solution]) + 1e-6
 end
-
-
-# seed = 0x946d4b7835e92ffa takes 90 minutes to solve!

src/MOI_bounded_oracle.jl

@@ -633,6 +633,8 @@ function solve(
     start_solution=nothing,
     fw_verbose=false,
     use_shadow_set=true,
+    custom_heuristics=[Heuristic()],
+    rounding_prob=1.0, 
     kwargs...,
 )
     blmo = convert(MathOptBLMO, lmo)
@@ -666,6 +668,8 @@ function solve(
         start_solution=start_solution,
         fw_verbose=fw_verbose,
         use_shadow_set=use_shadow_set,
+        custom_heuristics=custom_heuristics,
+        rounding_prob=rounding_prob,
         kwargs...,
     )
 end

src/callbacks.jl

@@ -86,32 +86,6 @@ function build_FW_callback(
             end
         end
 
-        if check_rounding_value && state.tt == FrankWolfe.pp
-            # round values
-            x_rounded = copy(state.x)
-            for idx in tree.branching_indices
-                x_rounded[idx] = round(state.x[idx])
-            end
-            # check linear feasibility
-            if is_linear_feasible(tree.root.problem.tlmo, x_rounded) &&
-               is_integer_feasible(tree, x_rounded)
-                # evaluate f(rounded)
-                val = tree.root.problem.f(x_rounded)
-                if val < tree.incumbent
-                    tree.root.updated_incumbent[] = true
-                    node = tree.nodes[tree.root.current_node_id[]]
-                    sol = FrankWolfeSolution(val, x_rounded, node, :rounded)
-                    push!(tree.solutions, sol)
-                    if tree.incumbent_solution === nothing ||
-                       sol.objective < tree.incumbent_solution.objective
-                        tree.incumbent_solution = sol
-                    end
-                    tree.incumbent = val
-                    Bonobo.bound!(tree, node.id)
-                end
-            end
-        end
-
         # check for time limit
         if isfinite(time_limit) && Dates.now() >= time_ref + Dates.Second(time_limit)
             return false

src/heuristics.jl

@@ -2,3 +2,180 @@
 
 ## Have a general interface for heuristics 
 ## so that the user can link in custom heuristics.
+"""
+    Boscia Heuristic
+
+Interface for heuristics in Boscia.    
+`h` is the heuristic function receiving as input the tree, the bounded LMO and a point x (the current node solution).
+It returns the heuristic solution (can be nothing, we check for that) and whether feasibility still has to be check.
+`prob` is the probability with which it will be called.        
+"""
+# Would 'Heuristic' also suffice? Or might we run into Identifer conflicts with other packages?
+struct Heuristic{F<:Function}
+    run_heuristic::F
+    prob::Float64
+    identifer::Symbol
+end
+
+# Create default heuristic. Doesn't do anything and should be called.
+Heuristic() = Heuristic(x -> nothing, -0.1, :default)
+
+"""
+Flip coin.
+"""
+function flip_coin(w=0.5, rng=Random.GLOBAL_RNG)
+    return rand(rng) ≤ w
+end
+
+"""
+Add a new solution found from the heuristic to the tree.
+"""
+function add_heuristic_solution(tree, x, val, heuristic_name::Symbol)
+    tree.root.updated_incumbent[] = true
+    node = tree.nodes[tree.root.current_node_id[]]
+    sol = FrankWolfeSolution(val, x, node, heuristic_name)
+    push!(tree.solutions, sol)
+    if tree.incumbent_solution === nothing ||
+        sol.objective < tree.incumbent_solution.objective
+        tree.incumbent_solution = sol
+    end
+    tree.incumbent = val
+    Bonobo.bound!(tree, node.id)
+end
+
+# TO DO: We might want to change the probability depending on the depth of the tree
+# or have other additional criteria on whether to run a heuristic
+"""
+Choose which heuristics to run by rolling a dice.
+"""
+function run_heuristics(tree, x, heuristic_list)
+    inner_lmo = tree.root.problem.tlmo.blmo
+    heuristic_lmo = TimeTrackingLMO(inner_lmo, tree.root.problem.integer_variables)
+
+    for heuristic in heuristic_list
+        if flip_coin(heuristic.prob)
+            list_x_heu, check_feasibility = heuristic.run_heuristic(tree, heuristic_lmo, x)
+
+            # check feasibility
+            if !isempty(list_x_heu)
+                min_val = Inf
+                min_idx = -1
+                for (i, x_heu) in enumerate(list_x_heu)
+                    feasible = check_feasibility ? is_linear_feasible(tree.root.problem.tlmo, x_heu) && is_integer_feasible(tree, x_heu) : true
+                    if feasible
+                        val = tree.root.problem.f(x_heu)
+                        if val < min_val
+                            min_val = val
+                            min_idx = i 
+                        end
+                    end
+                end
+
+                if min_val < tree.incumbent # Inf < Inf = false
+                    add_heuristic_solution(tree, list_x_heu[min_idx],min_val, heuristic.identifer)
+                end
+            end
+        end
+    end
+
+    # collect statistics from heuristic lmo
+    tree.root.options[:heu_ncalls] += heuristic_lmo.ncalls
+    return true
+end
+
+"""
+Simple rounding heuristic.
+"""
+function rounding_heuristic(tree::Bonobo.BnBTree, blmo::Boscia.TimeTrackingLMO, x)
+    x_rounded = copy(x)
+    for idx in tree.branching_indices
+        x_rounded[idx] = round(x[idx])
+    end
+    return [x_rounded], true
+end
+
+
+"""
+    follow-the-gradient
+Follow the gradient for a fixed number of steps and collect solutions on the way
+"""
+function follow_gradient_heuristic(tree::Bonobo.BnBTree, tlmo::Boscia.TimeTrackingLMO, x, k)
+    nabla = similar(x)
+    x_new = copy(x)
+    sols = []
+    for i in 1:k
+        tree.root.problem.g(nabla,x_new)
+        x_new = Boscia.compute_extreme_point(tlmo, nabla)
+        push!(sols, x_new)
+    end
+    return sols, false
+end
+
+
+"""
+Advanced lmo-aware rounding for binary vars. Rounding respecting the hidden feasible region structure.
+"""
+function rounding_lmo_01_heuristic(tree::Bonobo.BnBTree, tlmo::Boscia.TimeTrackingLMO, x)
+    nabla = zeros(length(x))
+    for idx in tree.branching_indices
+        nabla[idx] = 1 - 2*round(x[idx]) # (0.7, 0.3) -> (1, 0) -> (-1, 1) -> min -> (1,0)
+    end
+    x_rounded = Boscia.compute_extreme_point(tlmo, nabla)
+    return [x_rounded], false
+end
+
+"""
+Probability rounding for 0/1 problems.
+It decides based on the fractional value whether to ceil or floor the variable value. 
+Afterward, one call to Frank-Wolfe is performed to optimize the continuous variables.    
+"""
+function probability_rounding(tree::Bonobo.BnBTree, tlmo::Boscia.TimeTrackingLMO, x)
+    # save original bounds
+    node = tree.nodes[tree.root.current_node_id[]]
+    original_bounds = copy(node.local_bounds)
+
+    bounds = IntegerBounds()
+    for (i,x_i) in zip(tlmo.blmo.int_vars, x[tlmo.blmo.int_vars])
+        x_rounded = flip_coin(x_i) ? ceil(x_i) : floor(x_i)
+        push!(bounds, (i, x_rounded), :lessthan)
+        push!(bounds, (i, x_rounded), :greaterthan)
+    end
+
+    build_LMO(tlmo, tree.root.problem.integer_variable_bounds, bounds, tlmo.blmo.int_vars)
+
+    # check for feasibility and boundedness
+    status = check_feasibility(tlmo)
+    if status == INFEASIBLE || status == UNBOUNDED
+        @debug "LMO state in the probability rounding heuristic: $(status)"
+        # reset LMO to node state
+        build_LMO(tlmo, tree.root.problem.integer_variable_bounds, original_bounds, tlmo.blmo.int_vars)
+        # just return the point
+        return [x], false
+    end
+
+    v = compute_extreme_point(tlmo, rand(length(x)))
+    active_set = FrankWolfe.ActiveSet([(1.0, v)])
+
+    x_rounded, _, _, _ = solve_frank_wolfe(
+        tree.root.options[:variant],
+        tree.root.problem.f,
+        tree.root.problem.g,
+        tree.root.problem.tlmo,
+        active_set;
+        epsilon=node.fw_dual_gap_limit,
+        max_iteration=tree.root.options[:max_fw_iter],
+        line_search=tree.root.options[:lineSearch],
+        lazy=tree.root.options[:lazy],
+        lazy_tolerance=tree.root.options[:lazy_tolerance],
+        add_dropped_vertices=tree.root.options[:use_shadow_set],
+        use_extra_vertex_storage=tree.root.options[:use_shadow_set],
+        extra_vertex_storage=node.discarded_vertices,
+        callback=tree.root.options[:callback],
+        verbose=tree.root.options[:fwVerbose],
+    )
+
+    # reset LMO to node state
+    build_LMO(tlmo, tree.root.problem.integer_variable_bounds, original_bounds, tlmo.blmo.int_vars)
+
+    return [x_rounded], true
+end

src/interface.jl

@@ -51,6 +51,9 @@ fw_verbose             - if true, FrankWolfe logs are printed
 use_shadow_set         - The shadow set is the set of discarded vertices which is inherited by the children nodes.
                         It is used to avoid recomputing of vertices in case the LMO is expensive. In case of a cheap LMO,
                         performance might improve by disabling this option. 
+custom_heuristics      - List of  custom heuristic from the user.    
+prob_rounding          - The probability for calling the rounding heuristics. Since the feasibility has to be checked, it might
+                        expensive to do this for every node.                    
 """
 function solve(
     f,
@@ -84,6 +87,8 @@ function solve(
     start_solution=nothing,
     fw_verbose=false,
     use_shadow_set=true,
+    custom_heuristics=[Heuristic()],
+    rounding_prob=1.0,
     kwargs...,
 )
     if verbose
@@ -155,6 +160,9 @@ function solve(
         0.0,
     )
 
+    # Create standard heuristics
+    heuristics = vcat([Heuristic(rounding_heuristic, rounding_prob, :rounding)], custom_heuristics)
+
     Node = typeof(nodeEx)
     Value = Vector{Float64}
     tree = Bonobo.initialize(;
@@ -189,6 +197,8 @@ function solve(
                 :usePostsolve => use_postsolve,
                 :variant => variant,
                 :use_shadow_set => use_shadow_set,
+                :heuristics => heuristics, 
+                :heu_ncalls => 0,
             ),
         ),
         branch_strategy=branching_strategy,
@@ -530,6 +540,7 @@ function build_bnb_callback(
 
             result[:number_nodes] = tree.num_nodes
             result[:lmo_calls] = tree.root.problem.tlmo.ncalls
+            result[:heu_lmo_calls] = tree.root.options[:heu_ncalls]
             result[:list_num_nodes] = list_num_nodes_cb
             result[:list_lmo_calls_acc] = list_lmo_calls_cb
             result[:list_active_set_size] = list_active_set_size_cb
@@ -649,7 +660,13 @@ function postsolve(tree, result, time_ref, verbose, max_iteration_post)
         println("\t Dual Gap (relative): $(relative_gap(primal,tree_lb(tree)))\n")
         println("Search Statistics.")
         println("\t Total number of nodes processed: ", tree.num_nodes)
-        println("\t Total number of lmo calls: ", tree.root.problem.tlmo.ncalls)
+        if tree.root.options[:heu_ncalls] != 0
+            println("\t LMO calls over all nodes: ", tree.root.problem.tlmo.ncalls)
+            println("\t LMO calls in the heuristics: ", tree.root.options[:heu_ncalls])
+            println("\t Total number of lmo calls: ", tree.root.problem.tlmo.ncalls + tree.root.options[:heu_ncalls])
+        else
+            println("\t Total number of lmo calls: ", tree.root.problem.tlmo.ncalls)
+        end
         println("\t Total time (s): ", total_time_in_sec)
         println("\t LMO calls / sec: ", tree.root.problem.tlmo.ncalls / total_time_in_sec)
         println("\t Nodes / sec: ", tree.num_nodes / total_time_in_sec)