Heuristics Interface

src/MOI_bounded_oracle.jl

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

src/callbacks.jl

@@ -82,32 +82,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
-
         return true
     end
 end

src/heuristics.jl

@@ -2,3 +2,78 @@
 
 ## 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
+    h::Function
+    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)
+    return rand() ≤ w
+end
+
+"""
+Add a new solution found from the heuristic to the tree.
+"""
+function add_heuristic_solution(tree, x, val, heu::Symbol)
+    tree.root.updated_incumbent[] = true
+    node = tree.nodes[tree.root.current_node_id[]]
+    sol = FrankWolfeSolution(val, x, node, heu)
+    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
+
+"""
+Choose which heuristics to run by rolling a dice.
+"""
+# 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
+function run_heuristics(tree, x, heuristic_list)
+    for heuristic in heuristic_list
+        if flip_coin(heuristic.prob)
+            x_heu, check_feasibility = heuristic.h(tree, tree.root.problem.tlmo.blmo, x)
+
+            # check feasibility
+            if x_heu !== nothing 
+                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 < tree.incumbent
+                        add_heuristic_solution(tree, x_heu, val, heuristic.identifer)
+                    end
+                end
+            end
+        end
+    end
+end
+
+"""
+Simple rounding heuristic.
+"""
+function rounding_heuristic(tree::Bonobo.BnBTree, blmo::BoundedLinearMinimizationOracle, x)
+    x_rounded = copy(x)
+    for idx in tree.branching_indices
+        x_rounded[idx] = round(x[idx])
+    end
+    return x_rounded, true
+end

src/interface.jl

@@ -51,6 +51,7 @@ 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. The default is a dummy heuristic that return nothing and is, in practice, never called.                          
 """
 function solve(
     f,
@@ -84,6 +85,7 @@ function solve(
     start_solution=nothing,
     fw_verbose=false,
     use_shadow_set=true,
+    custom_heuristics=[Heuristic()],
     kwargs...,
 )
     if verbose
@@ -155,6 +157,9 @@ function solve(
         0.0,
     )
 
+    # Create standard heuristics
+    heuristics = vcat([Heuristic(rounding_heuristic, 0.7,:rounding)], custom_heuristics)
+
     Node = typeof(nodeEx)
     Value = Vector{Float64}
     tree = Bonobo.initialize(;
@@ -189,6 +194,7 @@ function solve(
                 :usePostsolve => use_postsolve,
                 :variant => variant,
                 :use_shadow_set => use_shadow_set,
+                :heuristics => heuristics, # Is a vector/list sufficent or would a dictonary be better?
             ),
         ),
         branch_strategy=branching_strategy,

src/managed_blmo.jl

@@ -260,6 +260,7 @@ function solve(
     start_solution=nothing,
     fw_verbose=false,
     use_shadow_set=true,
+    custom_heuristics=[Heuristic()],
     kwargs...,
 )
     blmo = ManagedBoundedLMO(sblmo, lower_bounds, upper_bounds, int_vars, n)
@@ -295,6 +296,7 @@ function solve(
         start_solution=start_solution,
         fw_verbose=fw_verbose,
         use_shadow_set=use_shadow_set,
+        custom_heuristics=custom_heuristics,
         kwargs...,
     )
 end

src/node.jl

@@ -303,6 +303,9 @@ function Bonobo.evaluate_node!(tree::Bonobo.BnBTree, node::FrankWolfeNode)
         @assert lower_bound <= tree.incumbent + 1e-5 "lower_bound <= tree.incumbent + 1e-5 : $(lower_bound) <= $(tree.incumbent)"
     end
 
+    # Call heuristic 
+    run_heuristics(tree, x, tree.root.options[:heuristics])
+
     return lower_bound, NaN
 end
 

src/polytope_blmos.jl

@@ -76,6 +76,56 @@
     return isapprox(sum(v), sblmo.N, atol=1e-4, rtol=1e-2)
 end
 
+"""
+Hyperplane aware rounding for the Probability simplex.
+"""
+function rounding_hyperplane_heuristic(tree::Bonobo.BnBTree, blmo::ManagedBoundedLMO{ProbabilitySimplexSimpleBLMO}, x) 
+    z = copy(x)
+    for idx in tree.branching_indices
+        z[idx] = round(x[idx])
+    end
+
+    N = blmo.simple_lmo.N
+    if sum(z) < N
+        while sum(z) < N
+            z = add_to_min(z, blmo.upper_bounds, tree.branching_indices)
+        end
+    elseif sum(z) > N
+        while sum(z) > N
+            z = remove_from_max(z, blmo.lower_bounds, tree.branching_indices)
+        end
+    end
+    return z, true
+end
+function add_to_min(x, ub, int_vars)
+    perm = sortperm(x)
+    j = findfirst(x->x != 0, x[perm])
+
+    for i in intersect(j:length(x), int_vars)
+        if x[perm[i]] < ub[perm[i]]
+            x[perm[i]] += 1
+            break
+        else
+            continue
+        end
+    end
+    return x
+end
+function remove_from_max(x, lb, int_vars)
+    perm = sortperm(x, rev = true)
+    j = findlast(x->x != 0, x[perm])
+
+    for i in intersect(1:j, int_vars)
+        if x[perm[i]] > lb[perm[i]]
+            x[perm[i]] -= 1
+            break
+        else
+            continue
+        end
+    end
+    return x
+end
+
 """
     UnitSimplexSimpleBLMO(N)
 
@@ -118,3 +168,21 @@
     end
     return sum(v) ≤ sblmo.N + 1e-3
 end
+
+"""
+Hyperplane aware rounding for the Unit simplex.
+"""
+function rounding_hyperplane_heuristic(tree::Bonobo.BnBTree, blmo::ManagedBoundedLMO{UnitSimplexSimpleBLMO}, x) 
+    z = copy(x)
+    for idx in tree.branching_indices
+        z[idx] = round(x[idx])
+    end
+
+    N = blmo.simple_lmo.N
+    if sum(z) > N
+        while sum(z) > N
+            z = remove_from_max(z, blmo.lower_bounds, tree.branching_indices)
+        end
+    end
+    return z, true
+end

test/LMO_test.jl

@@ -166,3 +166,50 @@ diffi = x_sol + 0.3*rand([-1,1], n)
     @test sum(isapprox.(x, x_sol, atol=1e-6, rtol=1e-2)) == n
     @test isapprox(f(x), f(result[:raw_solution]), atol=1e-6, rtol=1e-3)
 end
+
+n = 20
+x_sol = rand(1:Int(floor(n/4)), n)
+N = sum(x_sol)
+dir = vcat(fill(1, Int(floor(n/2))), fill(-1, Int(floor(n/2))), fill(0, mod(n,2)))
+diffi = x_sol + 0.3 * dir
+
+@testset "Custom Heuristic - Probability Simplex" begin
+    function f(x)
+        return 0.5 * sum((x[i] - diffi[i])^2 for i in eachindex(x))
+    end
+    function grad!(storage, x)
+        @. storage = x - diffi
+    end
+
+    sblmo = Boscia.ProbabilitySimplexSimpleBLMO(N)
+    heu = Boscia.Heuristic(Boscia.rounding_hyperplane_heuristic, 0.8, :hyperplane_rounding)
+
+    x, _, result =
+        Boscia.solve(f, grad!, sblmo, fill(0.0, n), fill(1.0*N, n), collect(1:n), n, custom_heuristics=[heu])
+
+    @test sum(isapprox.(x, x_sol, atol=1e-6, rtol=1e-2)) == n
+    @test isapprox(f(x), f(result[:raw_solution]), atol=1e-6, rtol=1e-3)
+end
+
+n = 20
+x_sol = rand(1:Int(floor(n/4)), n)
+diffi = x_sol + 0.3*rand([-1,1], n)
+
+@testset "Custom Heuristic - Unit Simplex" begin
+    function f(x)
+        return 0.5 * sum((x[i] - diffi[i])^2 for i in eachindex(x))
+    end
+    function grad!(storage, x)
+        @. storage = x - diffi
+    end
+
+    N = sum(x_sol) + floor(n/2)
+    sblmo = Boscia.UnitSimplexSimpleBLMO(N)
+    heu = Boscia.Heuristic(Boscia.rounding_hyperplane_heuristic, 0.8, :hyperplane_rounding)
+
+    x, _, result =
+        Boscia.solve(f, grad!, sblmo, fill(0.0, n), fill(1.0*N, n), collect(1:n), n, custom_heuristics=[heu])
+
+    @test sum(isapprox.(x, x_sol, atol=1e-6, rtol=1e-2)) == n
+    @test isapprox(f(x), f(result[:raw_solution]), atol=1e-6, rtol=1e-3)
+end

test/interface_test.jl

@@ -117,6 +117,30 @@ end
     end
 end
 
+@testset "Start with Active Set" begin
+
+    function f(x)
+        return 0.5 * sum((x[i] - diffi[i])^2 for i in eachindex(x))
+    end
+    function grad!(storage, x)
+        @. storage = x - diffi
+    end
+
+    int_vars = collect(1:n)
+    lbs = zeros(n)
+    ubs = ones(n)
+
+    sblmo = Boscia.CubeSimpleBLMO(lbs, ubs, int_vars)
+    direction =rand(n)
+    v = Boscia.bounded_compute_extreme_point(sblmo, direction, lbs, ubs, int_vars)
+    active_set = FrankWolfe.ActiveSet([(1.0, v)])
+
+    x, _, result = Boscia.solve(f, grad!, sblmo, lbs[int_vars], ubs[int_vars], int_vars, n, active_set=active_set)
+
+    @test x == round.(diffi)
+    @test isapprox(f(x), f(result[:raw_solution]), atol=1e-6, rtol=1e-3)
+end
+
 # Sparse Poisson regression
 # min_{w, b, z} ∑_i exp(w x_i + b) - y_i (w x_i + b) + α norm(w)^2
 # s.t. -N z_i <= w_i <= N z_i