Branching with non-trivial domain oracle (#207)

User has to provide a domain feasible point respecting the new bound constraints. We solve the projection problem.

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
 
 """
@@ -155,10 +153,13 @@ end
 """
 Find n linearly independent rows of A to build the starting point.
 """
-function linearly_independent_rows(A)
+function linearly_independent_rows(A; u=fill(1, size(A, 1)))
     S = []
     m, n = size(A)
     for i in 1:m
+        if iszero(u[i])
+            continue
+        end
         S_i= vcat(S, i)
         if rank(A[S_i,:])==length(S_i)
             S=S_i
@@ -178,8 +179,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 add_to_min function, we do not require x to have coordinates at zero.
+"""
+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 +208,6 @@ function remove_from_max(x)
         if x[perm[i]] > 1
             x[perm[i]] -= 1
             break
-        else
-            continue
         end
     end
     return x
@@ -237,6 +250,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, u=.!(iszero.(ub)))
+        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,23 +42,25 @@ 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)
@@ -72,6 +74,7 @@ verbose = true
         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,
     )

examples/strong_branching_portfolio.jl

@@ -2,7 +2,6 @@ using Boscia
 using FrankWolfe
 using Test
 using Random
-using SCIP
 using LinearAlgebra
 using Distributions
 import MathOptInterface
@@ -28,7 +27,7 @@ const Mi = (Ai + Ai') / 2
 @assert isposdef(Mi)
 
 function prepare_portfolio_lmo()
-    o = SCIP.Optimizer()
+    o = HiGHS.Optimizer()
     MOI.set(o, MOI.Silent(), true)
     MOI.empty!(o)
     x = MOI.add_variables(o, n)

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,