diff --git a/examples/nonlinear.jl b/examples/nonlinear.jl index 5d28281f1..4fac4e8f9 100644 --- a/examples/nonlinear.jl +++ b/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 diff --git a/examples/portfolio.jl b/examples/portfolio.jl index f8fdc8b34..8b0221384 100644 --- a/examples/portfolio.jl +++ b/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! diff --git a/src/MOI_bounded_oracle.jl b/src/MOI_bounded_oracle.jl index 3b7970a2b..f3e7f74fc 100644 --- a/src/MOI_bounded_oracle.jl +++ b/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 diff --git a/src/callbacks.jl b/src/callbacks.jl index 38b266346..dcb14f3ee 100644 --- a/src/callbacks.jl +++ b/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 diff --git a/src/heuristics.jl b/src/heuristics.jl index 16e672399..ae6d871e0 100644 --- a/src/heuristics.jl +++ b/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 diff --git a/src/interface.jl b/src/interface.jl index a1f36361c..2912eb12a 100644 --- a/src/interface.jl +++ b/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) diff --git a/src/managed_blmo.jl b/src/managed_blmo.jl index abeab4373..ccea6d51c 100644 --- a/src/managed_blmo.jl +++ b/src/managed_blmo.jl @@ -210,6 +210,19 @@ function build_LMO_correct(blmo::ManagedBoundedLMO, node_bounds) return true end +function check_feasibility(blmo::ManagedBoundedLMO) + for (lb,ub) in zip(blmo.lower_bounds, blmo.upper_bounds) + if ub < lb + return INFEASIBLE + end + end + return check_feasibility(blmo.simple_lmo, blmo.lower_bounds, blmo.upper_bounds, blmo.int_vars, blmo.n) +end + +function check_feasibility(simple_lmo::SimpleBoundableLMO, lb, ub, int_vars, n) + return true +end + # Check whether a split is valid, i.e. the upper and lower on variable vidx are not the same. function is_valid_split(tree::Bonobo.BnBTree, blmo::ManagedBoundedLMO, vidx::Int) idx = findfirst(x -> x == vidx, blmo.int_vars) @@ -260,6 +273,8 @@ function solve( start_solution=nothing, fw_verbose=false, use_shadow_set=true, + custom_heuristics=[Heuristic()], + rounding_prob=1.0, kwargs..., ) blmo = ManagedBoundedLMO(sblmo, lower_bounds, upper_bounds, int_vars, n) @@ -295,6 +310,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 diff --git a/src/node.jl b/src/node.jl index 9de63b69e..2fc5d9bbd 100644 --- a/src/node.jl +++ b/src/node.jl @@ -236,6 +236,7 @@ function Bonobo.evaluate_node!(tree::Bonobo.BnBTree, node::FrankWolfeNode) # check for feasibility and boundedness status = check_feasibility(tree.root.problem.tlmo) if status == INFEASIBLE + @debug "Problem at node $(node.id) infeasible" return NaN, NaN end if status == UNBOUNDED @@ -293,6 +294,7 @@ function Bonobo.evaluate_node!(tree::Bonobo.BnBTree, node::FrankWolfeNode) # Found an upper bound if is_integer_feasible(tree, x) node.ub = primal + @debug "Node $(node.id) has an integer solution." return lower_bound, primal # Sanity check: If the incumbent is better than the lower bound of the root node # and the root node is not integer feasible, something is off! @@ -302,6 +304,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 diff --git a/src/polytope_blmos.jl b/src/polytope_blmos.jl index d25ffece0..cd21416c0 100644 --- a/src/polytope_blmos.jl +++ b/src/polytope_blmos.jl @@ -62,7 +62,7 @@ function bounded_compute_extreme_point(sblmo::ProbabilitySimplexSimpleBLMO, d, l idx = findfirst(x -> x == i, int_vars) v[i] += min(ub[idx]-lb[idx], sblmo.N - sum(v)) else - v[i] += N - sum(v) + v[i] += sblmo.N - sum(v) end end return v @@ -76,6 +76,65 @@ function is_simple_linear_feasible(sblmo::ProbabilitySimplexSimpleBLMO, v) return isapprox(sum(v), sblmo.N, atol=1e-4, rtol=1e-2) end +function check_feasibility(sblmo::ProbabilitySimplexSimpleBLMO, lb, ub, int_vars, n) + m = n - length(int_vars) + if sum(lb) ≤ sblmo.N ≤ sum(ub) + m*sblmo.N + return OPTIMAL + else + INFEASIBLE + end +end + +""" +Hyperplane-aware rounding for the probability simplex. +""" +function rounding_hyperplane_heuristic(tree::Bonobo.BnBTree, tlmo::TimeTrackingLMO{ManagedBoundedLMO{ProbabilitySimplexSimpleBLMO}}, x) + z = copy(x) + for idx in tree.branching_indices + z[idx] = round(x[idx]) + end + + N = tlmo.blmo.simple_lmo.N + if sum(z) < N + while sum(z) < N + z = add_to_min(z, tlmo.blmo.upper_bounds, tree.branching_indices) + end + elseif sum(z) > N + while sum(z) > N + z = remove_from_max(z, tlmo.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 +177,29 @@ function is_simple_linear_feasible(sblmo::UnitSimplexSimpleBLMO, v) end return sum(v) ≤ sblmo.N + 1e-3 end + +function check_feasibility(sblmo::UnitSimplexSimpleBLMO, lb, ub, int_vars, n) + if sum(lb) ≤ sblmo.N + return OPTIMAL + else + INFEASIBLE + end +end + +""" +Hyperplane-aware rounding for the unit simplex. +""" +function rounding_hyperplane_heuristic(tree::Bonobo.BnBTree, tlmo::TimeTrackingLMO{ManagedBoundedLMO{UnitSimplexSimpleBLMO}}, x) + z = copy(x) + for idx in tree.branching_indices + z[idx] = round(x[idx]) + end + + N = tlmo.blmo.simple_lmo.N + if sum(z) > N + while sum(z) > N + z = remove_from_max(z, tlmo.blmo.lower_bounds, tree.branching_indices) + end + end + return [z], true +end diff --git a/test/LMO_test.jl b/test/LMO_test.jl index e1901a696..bfb284af7 100644 --- a/test/LMO_test.jl +++ b/test/LMO_test.jl @@ -123,9 +123,9 @@ end end n = 20 -x_sol = rand(1:Int(floor(n/4)), n) +x_sol = rand(1:floor(Int, 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))) +dir = vcat(fill(1, floor(Int, n/2)), fill(-1, floor(Int, n/2)), fill(0, mod(n,2))) diffi = x_sol + 0.3 * dir @testset "Probability Simplex LMO" begin @@ -146,7 +146,7 @@ diffi = x_sol + 0.3 * dir end n = 20 -x_sol = rand(1:Int(floor(n/4)), n) +x_sol = rand(1:floor(Int, n/4), n) diffi = x_sol + 0.3*rand([-1,1], n) @testset "Unit Simplex LMO" begin diff --git a/test/heuristics.jl b/test/heuristics.jl new file mode 100644 index 000000000..6005832c3 --- /dev/null +++ b/test/heuristics.jl @@ -0,0 +1,170 @@ +using Test +using Boscia +using FrankWolfe +using Random +using SCIP +import MathOptInterface +import Bonobo +using HiGHS +using Printf +using Dates +const MOI = MathOptInterface +const MOIU = MOI.Utilities + +n = 20 +x_sol = rand(1:floor(Int, n/4), n) +N = sum(x_sol) +dir = vcat(fill(1, floor(Int, n/2)), fill(-1, floor(Int, n/2)), fill(0, mod(n,2))) +diffi = x_sol + 0.3 * dir + +@testset "Hyperplane Aware Rounding - 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:floor(Int, n/4), n) +diffi = x_sol + 0.3*rand([-1,1], n) + +@testset "Hyperplane Aware Rounding - 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 + +@testset "Following Gradient 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) + depth = 5 + heu = Boscia.Heuristic((tree, blmo, x) -> Boscia.follow_gradient_heuristic(tree,blmo,x, depth), 1.0, :follow_gradient) + + x_heu, _, result_heu = + Boscia.solve(f, grad!, sblmo, fill(0.0, n), fill(1.0*N, n), collect(1:n), n, custom_heuristics=[heu]) + + x, _, result = Boscia.solve(f, grad!, sblmo, fill(0.0, n), fill(1.0*N, n), collect(1:n), 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) + + @test sum(isapprox.(x_heu, x_sol, atol=1e-6, rtol=1e-2)) == n + @test isapprox(f(x_heu), f(result_heu[:raw_solution]), atol=1e-6, rtol=1e-3) + + @test result[:lmo_calls] == result_heu[:lmo_calls] + @test result_heu[:heu_lmo_calls] > 0 +end + +@testset "Rounding 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) + + x_always, _, result_always = + Boscia.solve(f, grad!, sblmo, fill(0.0, n), fill(1.0*N, n), collect(1:n), n) + + sblmo = Boscia.UnitSimplexSimpleBLMO(N) + x, _, result = Boscia.solve(f, grad!, sblmo, fill(0.0, n), fill(1.0*N, n), collect(1:n), n, rounding_prob = 0.5) + + @test sum(isapprox.(x_always, x_sol, atol=1e-6, rtol=1e-2)) == 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 = 30 +diffi = Random.rand(Bool, n) * 0.6 .+ 0.3 + +@testset "Probability Rounding - Unit Cube" 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 + + lbs = zeros(n) + ubs = ones(n) + int_vars = unique!(rand(1:n, floor(Int, n/2))) + x_sol = copy(diffi) + x_sol[int_vars] = round.(x_sol[int_vars]) + + sblmo = Boscia.CubeSimpleBLMO(lbs, ubs, int_vars) + heu = Boscia.Heuristic(Boscia.probability_rounding, 0.6, :probability_rounding) + + x, _, result = + Boscia.solve(f, grad!, sblmo, lbs[int_vars], ubs[int_vars], int_vars, n, custom_heuristics=[heu], rounding_prob=0.0) + + @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 = round.(rand(n)) +N = sum(x_sol) +dir = sign.(iszero.(x_sol) .- 0.5) +diffi = x_sol + 0.3 * dir + +@testset "Probability rounding - 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 + + int_vars = unique!(rand(1:n, floor(Int, n/2))) + m = length(int_vars) + cont_vars = setdiff(collect(1:n), int_vars) + x_sol[cont_vars] = diffi[cont_vars] + + sblmo = Boscia.ProbabilitySimplexSimpleBLMO(N) + heu = Boscia.Heuristic(Boscia.probability_rounding, 0.6, :probability_rounding) + + x, _, result = + Boscia.solve(f, grad!, sblmo, fill(0.0, m), fill(1.0, m), int_vars, n, custom_heuristics=[heu], rounding_prob=0.0) + + @test f(x) ≥ f(x_sol) + if isapprox(sum(x_sol), N) + @test isapprox(f(x), f(x_sol)) + end + @test sum(isapprox.(x[int_vars], x_sol[int_vars], atol=1e-6, rtol=1e-2)) == m + @test isapprox(f(x), f(result[:raw_solution]), atol=1e-6, rtol=1e-3) +end diff --git a/test/interface_test.jl b/test/interface_test.jl index f97fe14a3..572f0b557 100644 --- a/test/interface_test.jl +++ b/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 diff --git a/test/runtests.jl b/test/runtests.jl index 8df50d075..c9ed78042 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -35,6 +35,7 @@ end include("LMO_test.jl") include("indicator_test.jl") +include("heuristics.jl") # Takes pretty long, only include if you want to test this specifically