Branching with non-trivial domain oracle

examples/oed_utils.jl

@@ -6,10 +6,8 @@ m    - number of experiments.
 fusion - boolean deiciding whether we build the fusion or standard problem.
 corr - boolean deciding whether we build the independent or correlated data.   
 """
-function build_data(m)
+function build_data(m, n)
     # set up
-    n = Int(floor(m/10))
-
     B = rand(m,n)
     B = B'*B
     @assert isposdef(B)
@@ -22,7 +20,7 @@ function build_data(m)
     u = floor(N/3)
     ub = rand(1.0:u, m)
 
-    return A, n, N, ub 
+    return A, N, ub 
 end
 
 """
@@ -178,8 +176,22 @@ function add_to_min(x, u)
         if x[perm[i]] < u[perm[i]]
             x[perm[i]] += 1
             break
-        else
-            continue
+        end
+    end
+    return x
+end
+
+"""
+We want to add to the smallest value of x while respecting the upper bounds u.
+In constrast to the function above, we do not require x to have zero entries.
+"""
+function add_to_min2(x,u)
+    perm = sortperm(x)
+
+    for i in perm
+        if x[i] < u[i]
+            x[i] += 1
+            break
         end
     end
     return x
@@ -193,8 +205,6 @@ function remove_from_max(x)
         if x[perm[i]] > 1
             x[perm[i]] -= 1
             break
-        else
-            continue
         end
     end
     return x
@@ -237,6 +247,52 @@ function build_start_point(A, N, ub)
     return x, active_set, S
 end
 
+"""
+Build domain feasible point for any node.
+The point has to be domain feasible as well as respect the 
+node bounds. If no such point can be constructed, return nothing.
+"""
+function build_domain_point_function(domain_oracle, A, N, int_vars, initial_lb, initial_ub)
+    return function domain_point(local_bounds)
+        lb = initial_lb
+        ub = initial_ub
+        for idx in int_vars
+            if haskey(local_bounds.lower_bounds, idx)
+                lb[idx] = max(initial_lb[idx], local_bounds.lower_bounds[idx])
+            end
+            if haskey(local_bounds.upper_bounds, idx)
+                ub[idx] = min(initial_ub[idx], local_bounds.upper_bounds[idx])
+            end
+        end
+        # Node itself infeasible
+        if sum(lb) > N 
+            return nothing
+        end
+        # No intersection between node and domain
+        if !domain_oracle(ub)
+            return nothing
+        end
+        x = lb
+
+        S = linearly_independent_rows(A)
+        while sum(x) <= N
+            if sum(x) == N 
+                if domain_oracle(x)
+                    return x 
+                else 
+                    return nothing 
+                end 
+            end
+            if iszero(x[S]-ub[S])
+                y = add_to_min2(x[S], ub[S])
+                x[S] = y
+            else
+                x = add_to_min2(x, ub)
+            end
+        end
+    end
+end
+
 """
 Create first incumbent for Boscia and custom BB in a greedy fashion.
 """

examples/optimal_experiment_design.jl

@@ -42,36 +42,39 @@ include("oed_utils.jl")
 
 
 m = 50
+n = Int(floor(m / 10))
 verbose = true
 
 ## A-Optimal Design Problem
 @testset "A-Optimal Design" begin
 
-    Ex_mat, n, N, ub = build_data(m)
+    Ex_mat, N, ub = build_data(m, n)
 
     g, grad! = build_a_criterion(Ex_mat, build_safe=false)
     blmo = build_blmo(m, N, ub)
     heu = Boscia.Heuristic(Boscia.rounding_hyperplane_heuristic, 0.7, :hyperplane_aware_rounding)
     domain_oracle = build_domain_oracle(Ex_mat, n)
+    domain_point = build_domain_point_function(domain_oracle, Ex_mat, N, collect(1:m), fill(0.0, m), ub)
 
     # precompile
     line_search = FrankWolfe.MonotonicGenericStepsize(FrankWolfe.Adaptive(), domain_oracle)
     x0, active_set = build_start_point(Ex_mat, N, ub)
     z = greedy_incumbent(Ex_mat, N, ub)
-    x, _, _ = Boscia.solve(g, grad!, blmo, active_set=active_set, start_solution=z, time_limit=10, verbose=false, domain_oracle=domain_oracle, custom_heuristics=[heu], line_search=line_search)
+    x, _, _ = Boscia.solve(g, grad!, blmo, active_set=active_set, start_solution=z, time_limit=10, verbose=false, domain_oracle=domain_oracle, find_domain_point=domain_point, custom_heuristics=[heu], line_search=line_search)
 
     # proper run with MGLS and Adaptive
     line_search = FrankWolfe.MonotonicGenericStepsize(FrankWolfe.Adaptive(), domain_oracle)
     x0, active_set = build_start_point(Ex_mat, N, ub)
     z = greedy_incumbent(Ex_mat, N, ub)
-    x, _, result = Boscia.solve(
+   x, _, result = Boscia.solve(
         g, 
         grad!, 
         blmo, 
         active_set=active_set, 
         start_solution=z, 
         verbose=verbose, 
         domain_oracle=domain_oracle, 
+        find_domain_point=domain_point, 
         custom_heuristics=[heu], 
         line_search=line_search,
     )
@@ -89,29 +92,31 @@ verbose = true
         start_solution=z, 
         verbose=verbose, 
         domain_oracle=domain_oracle, 
+        find_domain_point=domain_point, 
         custom_heuristics=[heu], 
         line_search=line_search,
     )
 
-    @test result_s[:dual_bound] <= g(x) + 1e-4
-    @test result[:dual_bound] <= g(x_s) + 1e-4
+    @test result_s[:dual_bound] <= g(x) + 1e-3
+    @test result[:dual_bound] <= g(x_s) + 1e-3
     @test isapprox(g(x), g(x_s), atol=1e-3) 
 end 
 
 ## D-Optimal Design Problem
 @testset "D-optimal Design" begin
-    Ex_mat, n, N, ub = build_data(m)
+    Ex_mat, N, ub = build_data(m, n)
 
     g, grad! = build_d_criterion(Ex_mat, build_safe=false)
     blmo = build_blmo(m, N, ub)
     heu = Boscia.Heuristic(Boscia.rounding_hyperplane_heuristic, 0.7, :hyperplane_aware_rounding)
     domain_oracle = build_domain_oracle(Ex_mat, n)
+    domain_point = build_domain_point_function(domain_oracle, Ex_mat, N, collect(1:m), fill(0.0, m), ub)
 
     # precompile
     line_search = FrankWolfe.MonotonicGenericStepsize(FrankWolfe.Adaptive(), domain_oracle)
     x0, active_set = build_start_point(Ex_mat, N, ub)
     z = greedy_incumbent(Ex_mat, N, ub)
-    x, _, _ = Boscia.solve(g, grad!, blmo, active_set=active_set, start_solution=z, time_limit=10, verbose=false, domain_oracle=domain_oracle, custom_heuristics=[heu], line_search=line_search)
+    x, _, _ = Boscia.solve(g, grad!, blmo, active_set=active_set, start_solution=z, time_limit=10, verbose=false, domain_oracle=domain_oracle, find_domain_point=domain_point, custom_heuristics=[heu], line_search=line_search)
 
     # proper run with MGLS and Adaptive
     line_search = FrankWolfe.MonotonicGenericStepsize(FrankWolfe.Adaptive(), domain_oracle)
@@ -125,6 +130,7 @@ end
         start_solution=z, 
         verbose=verbose, 
         domain_oracle=domain_oracle, 
+        find_domain_point=domain_point, 
         custom_heuristics=[heu], 
         line_search=line_search,
     )
@@ -142,6 +148,7 @@ end
         start_solution=z, 
         verbose=verbose, 
         domain_oracle=domain_oracle, 
+        find_domain_point=domain_point, 
         custom_heuristics=[heu], 
         line_search=line_search,
     )

src/Boscia.jl

@@ -4,6 +4,7 @@ using FrankWolfe
 import FrankWolfe: compute_extreme_point
 export compute_extreme_point
 using Random
+using LinearAlgebra
 import Bonobo
 using Printf
 using Dates

src/MOI_bounded_oracle.jl

@@ -631,7 +631,8 @@ function solve(
     strong_convexity=0.0,
     sharpness_constant = 0.0,
     sharpness_exponent = Inf,
-    domain_oracle=x -> true,
+    domain_oracle=_trivial_domain,
+    find_domain_point= _trivial_domain_point,
     start_solution=nothing,
     fw_verbose=false,
     use_shadow_set=true,
@@ -671,6 +672,7 @@ function solve(
         sharpness_constant=sharpness_constant,
         sharpness_exponent=sharpness_exponent,
         domain_oracle=domain_oracle,
+        find_domain_point=find_domain_point,
         start_solution=start_solution,
         fw_verbose=fw_verbose,
         use_shadow_set=use_shadow_set,

src/heuristics.jl

@@ -53,7 +53,7 @@ function run_heuristics(tree, x, heuristic_list; rng=Random.GLOBAL_RNG)
                 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) : false
+                    feasible = check_feasibility ? is_linear_feasible(tree.root.problem.tlmo, x_heu) && is_integer_feasible(tree, x_heu) && tree.root.options[:domain_oracle](x_heu) : false
                     if feasible
                         val = tree.root.problem.f(x_heu)
                         if val < min_val

src/interface.jl

@@ -49,8 +49,11 @@ sharpness_constant     - the constant M > 0 for (θ, M)-sharpness.
                          where X* is the set minimizer of f. 
 sharpness_exponent     - the exponent θ ∈ [0, 1/2] for (θ, M)-sharpness.
 domain_oracle          - For a point x: returns true if x is in the domain of f, else false. Per default is true.
-                         In case of the non trivial domain oracle, the starting point has to be feasible for f. Also, depending 
+                         In case of the non trivial domain oracle, the starting point has to be feasible for f. Additionally,
+                         the user has to provide a function `domain_point`, see below. Also, depending 
                          on the Line Search method, you might have to provide the domain oracle to it, too.
+find_domain_point      - Given the current node bounds return a domain feasible point respecting the bounds.
+                         If no such point can be found, return nothing.                       
 start_solution         - initial solution to start with an incumbent
 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.
@@ -98,7 +101,8 @@ function solve(
     strong_convexity=0.0,
     sharpness_constant = 0.0,
     sharpness_exponent = Inf,
-    domain_oracle=x -> true,
+    domain_oracle=_trivial_domain,
+    find_domain_point= _trivial_domain_point,
     start_solution=nothing,
     fw_verbose=false,
     use_shadow_set=true,
@@ -146,6 +150,10 @@ function solve(
 
     global_bounds = build_global_bounds(blmo, integer_variables)
 
+    if typeof(domain_oracle) != typeof(_trivial_domain) && typeof(find_domain_point) == typeof(_trivial_domain_point)
+        @warn "For a non trivial domain oracle, please provide the DOMAIN POINT function. Otherwise, Boscia might not converge."
+    end
+
     v = []
     if active_set === nothing
         direction = collect(1.0:n)
@@ -201,6 +209,7 @@ function solve(
             global_tightenings=IntegerBounds(),
             options=Dict{Symbol,Any}(
                 :domain_oracle => domain_oracle,
+                :find_domain_point => find_domain_point,
                 :dual_gap => dual_gap,
                 :dual_gap_decay_factor => dual_gap_decay_factor,
                 :dual_tightening => dual_tightening,

src/managed_blmo.jl

@@ -271,7 +271,8 @@ function solve(
     strong_convexity=0.0,
     sharpness_constant = 0.0,
     sharpness_exponent = Inf,
-    domain_oracle=x -> true,
+    domain_oracle=_trivial_domain,
+    find_domain_point= _trivial_domain_point,
     start_solution=nothing,
     fw_verbose=false,
     use_shadow_set=true,
@@ -313,6 +314,7 @@ function solve(
         sharpness_constant=sharpness_constant,
         sharpness_exponent=sharpness_exponent,
         domain_oracle=domain_oracle,
+        find_domain_point=find_domain_point,
         start_solution=start_solution,
         fw_verbose=fw_verbose,
         use_shadow_set=use_shadow_set,

src/node.jl

@@ -85,10 +85,30 @@ function Bonobo.get_branching_nodes_info(tree::Bonobo.BnBTree, node::FrankWolfeN
         split_vertices_set!(node.discarded_vertices, tree, vidx, x, node.local_bounds)
 
     # Sanity check
-    @assert isapprox(sum(active_set_left.weights), 1.0)
+    @assert isapprox(sum(active_set_left.weights), 1.0) "sum weights left: $(sum(active_set_left.weights))"
     @assert sum(active_set_left.weights .< 0) == 0
-    @assert isapprox(sum(active_set_right.weights), 1.0)
+    for v in active_set_left.atoms
+        if !(v[vidx] <= floor(x[vidx]) + tree.options.atol)
+            error("active_set_left\n$(v)\n$vidx, $(x[vidx]), $(v[vidx])")
+        end
+    end
+    @assert isapprox(sum(active_set_right.weights), 1.0) "sum weights right: $(sum(active_set_right.weights))"
     @assert sum(active_set_right.weights .< 0) == 0
+    for v in active_set_right.atoms
+        if !(v[vidx] >= ceil(x[vidx]) - tree.options.atol)
+            error("active_set_right\n$(v)\n$vidx, $(x[vidx]), $(v[vidx])")
+        end
+    end
+    for v in discarded_set_left.storage
+        if !(v[vidx] <= floor(x[vidx]) + tree.options.atol)
+            error("storage left\n$(v)\n$vidx, $(x[vidx]), $(v[vidx])")
+        end
+    end
+    for v in discarded_set_right.storage
+        if !(v[vidx] >= ceil(x[vidx]) - tree.options.atol)
+            error("storage right\n$(v)\n$vidx, $(x[vidx]), $(v[vidx])")
+        end
+    end
 
     # add new bounds to the feasible region left and right
     # copy bounds from parent
@@ -108,26 +128,15 @@ function Bonobo.get_branching_nodes_info(tree::Bonobo.BnBTree, node::FrankWolfeN
     fw_dual_gap_limit = tree.root.options[:dual_gap_decay_factor] * node.fw_dual_gap_limit
     fw_dual_gap_limit = max(fw_dual_gap_limit, tree.root.options[:min_node_fw_epsilon])
 
-    # Sanity check
-    for v in active_set_left.atoms
-        if !(v[vidx] <= floor(x[vidx]) + tree.options.atol)
-            error("active_set_left\n$(v)\n$vidx, $(x[vidx]), $(v[vidx])")
-        end
-    end
-    for v in discarded_set_left.storage
-        if !(v[vidx] <= floor(x[vidx]) + tree.options.atol)
-            error("storage left\n$(v)\n$vidx, $(x[vidx]), $(v[vidx])")
-        end
-    end
-    for v in active_set_right.atoms
-        if !(v[vidx] >= ceil(x[vidx]) - tree.options.atol)
-            error("active_set_right\n$(v)\n$vidx, $(x[vidx]), $(v[vidx])")
-        end
+    # in case of non trivial domain oracle: Only split if the iterate is still domain feasible
+    x_left = FrankWolfe.compute_active_set_iterate!(active_set_left) 
+    x_right = FrankWolfe.compute_active_set_iterate!(active_set_right)
+
+    if !tree.root.options[:domain_oracle](x_left)
+        active_set_left = build_active_set_by_domain_oracle(active_set_left, tree, varbounds_left, node)
     end
-    for v in discarded_set_right.storage
-        if !(v[vidx] >= ceil(x[vidx]) - tree.options.atol)
-            error("storage right\n$(v)\n$vidx, $(x[vidx]), $(v[vidx])")
-        end
+    if !tree.root.options[:domain_oracle](x_right)
+        active_set_right = build_active_set_by_domain_oracle(active_set_right, tree, varbounds_right, node)
     end
 
     # update the LMO
@@ -156,13 +165,8 @@ function Bonobo.get_branching_nodes_info(tree::Bonobo.BnBTree, node::FrankWolfeN
         dual_gap=NaN,
     )
 
-    # in case of non trivial domain oracle: Only split if the iterate is still domain feasible
-    x_left = FrankWolfe.compute_active_set_iterate!(active_set_left)
-    x_right = FrankWolfe.compute_active_set_iterate!(active_set_right)
-    domain_oracle = tree.root.options[:domain_oracle]
-
-    domain_right = domain_oracle(x_right)
-    domain_left = domain_oracle(x_left)
+    domain_right = !isempty(active_set_right)
+    domain_left = !isempty(active_set_left)
 
     nodes = if !prune_left && !prune_right && domain_right && domain_left
         [node_info_left, node_info_right]

src/problem.jl

@@ -8,6 +8,16 @@ Enum for the solving stage
     TIME_LIMIT_REACHED = 3
 end
 
+"""
+Trivial domain function.
+"""
+_trivial_domain(x) = true
+
+"""
+Trivial domain point function.
+"""
+_trivial_domain_point(local_bounds::IntegerBounds) = nothing
+
 abstract type AbstractSimpleOptimizationProblem end
 
 """