From 632ca6c1d126f8c22e5b0d56b691829e597330b5 Mon Sep 17 00:00:00 2001
From: dhendryc <92737336+dhendryc@users.noreply.github.com>
Date: Fri, 8 Nov 2024 13:23:34 +0100
Subject: [PATCH] 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                  | 77 ++++++++++++++++++---
 examples/optimal_experiment_design.jl  | 19 +++--
 examples/strong_branching_portfolio.jl |  3 +-
 src/Boscia.jl                          |  1 +
 src/MOI_bounded_oracle.jl              |  4 +-
 src/heuristics.jl                      |  2 +-
 src/interface.jl                       | 13 +++-
 src/managed_blmo.jl                    |  4 +-
 src/node.jl                            | 60 ++++++++--------
 src/problem.jl                         | 10 +++
 src/utilities.jl                       | 96 ++++++++++++++++++++++++++
 test/runtests.jl                       |  1 +
 12 files changed, 240 insertions(+), 50 deletions(-)

diff --git a/examples/oed_utils.jl b/examples/oed_utils.jl
index 32064d071..ce8d9495c 100644
--- a/examples/oed_utils.jl
+++ b/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.
 """
diff --git a/examples/optimal_experiment_design.jl b/examples/optimal_experiment_design.jl
index d9078d0c5..9973966ce 100644
--- a/examples/optimal_experiment_design.jl
+++ b/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,
     )
diff --git a/examples/strong_branching_portfolio.jl b/examples/strong_branching_portfolio.jl
index 577848163..ab01912fe 100644
--- a/examples/strong_branching_portfolio.jl
+++ b/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)
diff --git a/src/Boscia.jl b/src/Boscia.jl
index 6853fa541..6bb6a11c5 100644
--- a/src/Boscia.jl
+++ b/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
diff --git a/src/MOI_bounded_oracle.jl b/src/MOI_bounded_oracle.jl
index 39e883401..89b3f20d2 100644
--- a/src/MOI_bounded_oracle.jl
+++ b/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,
diff --git a/src/heuristics.jl b/src/heuristics.jl
index c30392a8c..0a0e4d5ab 100644
--- a/src/heuristics.jl
+++ b/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
diff --git a/src/interface.jl b/src/interface.jl
index 772d9c9a0..e08054c92 100644
--- a/src/interface.jl
+++ b/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,
diff --git a/src/managed_blmo.jl b/src/managed_blmo.jl
index c64967e28..3a580c716 100644
--- a/src/managed_blmo.jl
+++ b/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,
diff --git a/src/node.jl b/src/node.jl
index 1a201904e..1de2146fb 100644
--- a/src/node.jl
+++ b/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]
diff --git a/src/problem.jl b/src/problem.jl
index 8d1479c72..10fd0cba9 100644
--- a/src/problem.jl
+++ b/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
 
 """
diff --git a/src/utilities.jl b/src/utilities.jl
index d56efdb13..afc7dafab 100644
--- a/src/utilities.jl
+++ b/src/utilities.jl
@@ -155,6 +155,102 @@ function split_vertices_set!(
     return (discarded_set, right_as)
 end
 
+"""
+Build a new start point and active set in case the split active set
+does not lead to a domain feasible iterate.
+First, try filtering the active set by the domain oracle.
+If all vertices are domain infeasible, solve the projection problem
+1/2 * ||x - x*||_2^2 
+where x* is a domain- and bound-feasible point provided by the user.
+"""
+function build_active_set_by_domain_oracle(
+    active_set::FrankWolfe.ActiveSet{T,R},
+    tree,
+    local_bounds::IntegerBounds,
+    node;
+    atol=1e-5,
+    rtol=1e-5,
+) where {T,R}
+    # Check if node problem is even feasible
+    build_LMO(
+        tree.root.problem.tlmo,
+        tree.root.problem.integer_variable_bounds,
+        local_bounds,
+        tree.root.problem.integer_variables,
+    )
+    status = check_feasibility(tree.root.problem.tlmo)
+    if status == INFEASIBLE
+        build_LMO(
+            tree.root.problem.tlmo,
+            tree.root.problem.integer_variable_bounds,
+            node.local_bounds,
+            tree.root.problem.integer_variables,
+        )
+        active_set.empty!()
+        return active_set
+    end
+    # Filtering
+    del_indices = BitSet()
+    for (idx, tup) in enumerate(active_set)
+        (λ, a) = tup
+        if !tree.root.options[:domain_oracle](a)
+            push!(del_indices, idx)
+        end
+    end
+    # At least one vertex is domain feasible.
+    if length(del_indices) < length(active_set.weights)
+        deleteat!(active_set, del_indices)
+        @assert !isempty(active_set)
+        FrankWolfe.active_set_renormalize!(active_set)
+        FrankWolfe.compute_active_set_iterate!(active_set)
+    # No vertex is domain feasible
+    else
+        x_star = tree.root.options[:find_domain_point](local_bounds)
+        # No domain feasible point can be build.
+        # Node can be pruned.
+        if x_star === nothing
+            deleteat!(active_set, del_indices)
+        else
+            inner_f(x) = 1/2 * LinearAlgebra.norm(x - x_star)^2
+
+            function inner_grad!(storage, x)
+                storage .= x - x_star
+                return storage
+            end
+
+            function build_inner_callback(tree)
+                return function inner_callback(state, active_set, kwargs...)
+                    # stop as soon as we find a domain feasible point.
+                    if tree.root.options[:domain_oracle](state.x)
+                        return false
+                    end
+                end
+            end
+            inner_callback = build_inner_callback(tree)
+
+            x, _, _, _, _, active_set = FrankWolfe.blended_pairwise_conditional_gradient(
+                inner_f,
+                inner_grad!,
+                tree.root.problem.tlmo,
+                active_set,
+                callback=inner_callback,
+                lazy=true,
+            )
+            @assert tree.root.options[:domain_oracle](x)
+        end
+    end
+    build_LMO(
+        tree.root.problem.tlmo,
+        tree.root.problem.integer_variable_bounds,
+        node.local_bounds,
+        tree.root.problem.integer_variables,
+    )
+    return active_set
+end
+
+"""
+Check if a given point v satisfies the given bounds.
+"""
 function is_bound_feasible(bounds::IntegerBounds, v; atol=1e-5)
     for (idx, set) in bounds.lower_bounds
         if v[idx] < set - atol
diff --git a/test/runtests.jl b/test/runtests.jl
index c9ed78042..6ac70b25c 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -45,6 +45,7 @@ include("sparse_regression.jl")
 include("poisson.jl")
 include("mean_risk.jl")
 include("time_limit.jl")
+include("strong_convexity_and_sharpness.jl")
 
 n = 10
 const diff1 = rand(Bool, n) * 0.8 .+ 1.1