Add BDICG and rewrite the generation of the starting point in DICG/BDICG

examples/oed_utils.jl

@@ -259,8 +259,8 @@ 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
+        lb = copy(initial_lb)
+        ub = copy(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])

src/frank_wolfe_variants.jl

@@ -10,6 +10,7 @@ It may also implement `build_frank_wolfe_workspace(x)` which creates a
 workspace structure that is passed as last argument to `solve_frank_wolfe`.
 """
 abstract type FrankWolfeVariant end
+abstract type DecompositionInvariant <: FrankWolfeVariant end
 
 # default printing for FrankWolfeVariant is just showing the type
 Base.print(io::IO, fw::FrankWolfeVariant) = print(io, split(string(typeof(fw)), ".")[end])
@@ -185,24 +186,21 @@ Base.print(io::IO, ::BlendedPairwiseConditionalGradient) =
 The Decomposition-invariant Frank-Wolfe. 
 
 """
-struct DecompositionInvariantConditionalGradient <: FrankWolfeVariant
+struct DecompositionInvariantConditionalGradient <: DecompositionInvariant
     use_strong_lazy::Bool
     use_DICG_warm_start::Bool
     use_strong_warm_start::Bool
-    build_dicg_start_point::Function
 end
 
 function DecompositionInvariantConditionalGradient(;
     use_strong_lazy=false,
     use_DICG_warm_start=false,
     use_strong_warm_start=false,
-    build_dicg_start_point=trivial_build_dicg_start_point,
 )
     return DecompositionInvariantConditionalGradient(
         use_strong_lazy,
         use_DICG_warm_start,
         use_strong_warm_start,
-        build_dicg_start_point,
     )
 end
 
@@ -228,32 +226,7 @@ function solve_frank_wolfe(
     domain_oracle=_trivial_domain,
     kwargs...,
 )
-    # We keep track of computed extreme points by creating logging callback.
-    function make_callback(pre_computed_set)
-        return function DICG_callback(state, kwargs...)
-            if !callback(state, pre_computed_set)
-                return false
-            end
-            return true
-        end
-    end
-
-    x0 = dicg_start_point_initialize(
-        lmo,
-        active_set,
-        pre_computed_set,
-        frank_wolfe_variant.build_dicg_start_point;
-        domain_oracle=domain_oracle,
-    )
-
-    if x0 === nothing || !domain_oracle(x0)
-        return NaN, Inf, Inf, pre_computed_set
-    else
-        @assert is_linear_feasible(lmo, x0)
-    end
-
-    # In case of the postprocessing, no callback is provided.
-    DICG_callback = callback !== nothing ? make_callback(pre_computed_set) : nothing
+    x0, DICG_callback = init_decomposition_invariant_state(active_set, pre_computed_set, callback)
 
     x, _, primal, dual_gap, status, _ = FrankWolfe.decomposition_invariant_conditional_gradient(
         f,
@@ -272,25 +245,91 @@ function solve_frank_wolfe(
         callback=DICG_callback,
         extra_vertex_storage=pre_computed_set,
     )
-    if pre_computed_set !== nothing
-        if frank_wolfe_variant.use_strong_warm_start
-            indices_to_delete = []
-            for idx in eachindex(pre_computed_set)
-                atom = pre_computed_set[idx]
-                if !is_inface_feasible(lmo, atom, x)
-                    push!(indices_to_delete, idx)
-                end
-            end
-            deleteat!(pre_computed_set, indices_to_delete)
-        end
-    end
+
+    clean_up_pre_computed_set!(lmo, pre_computed_set, x, frank_wolfe_variant)
 
     return x, primal, dual_gap, status, pre_computed_set
 end
 
 Base.print(io::IO, ::DecompositionInvariantConditionalGradient) =
     print(io, "Decompostion-Invariant-Frank-Wolfe")
 
+"""
+   BDICG-Frank-Wolfe
+
+The Blended Decomposition-invariant Frank-Wolfe. 
+
+"""
+struct BlendedDecompositionInvariantConditionalGradient <: DecompositionInvariant
+    use_strong_lazy::Bool
+    use_BDICG_warm_start::Bool
+    use_strong_warm_start::Bool
+end
+
+function BlendedDecompositionInvariantConditionalGradient(;
+    use_strong_lazy=false,
+    use_BDICG_warm_start=false,
+    use_strong_warm_start=false,
+)
+    return BlendedDecompositionInvariantConditionalGradient(
+        use_strong_lazy,
+        use_BDICG_warm_start,
+        use_strong_warm_start,
+    )
+end
+
+function solve_frank_wolfe(
+    frank_wolfe_variant::BlendedDecompositionInvariantConditionalGradient,
+    f,
+    grad!,
+    lmo,
+    active_set;
+    line_search::FrankWolfe.LineSearchMethod=FrankWolfe.Secant(),
+    epsilon=1e-7,
+    max_iteration=10000,
+    add_dropped_vertices=false,
+    use_extra_vertex_storage=false,
+    extra_vertex_storage=nothing,
+    callback=nothing,
+    lazy=false,
+    lazy_tolerance=2.0,
+    timeout=Inf,
+    verbose=false,
+    workspace=nothing,
+    pre_computed_set=nothing,
+    domain_oracle=_trivial_domain,
+    kwargs...,
+)
+
+    x0, BDICG_callback = init_decomposition_invariant_state(active_set, pre_computed_set, callback)
+
+    x, _, primal, dual_gap, status, _ =
+        FrankWolfe.blended_decomposition_invariant_conditional_gradient(
+            f,
+            grad!,
+            lmo,
+            x0;
+            line_search=line_search,
+            epsilon=epsilon,
+            max_iteration=max_iteration,
+            verbose=verbose,
+            timeout=timeout,
+            lazy=lazy,
+            use_strong_lazy=frank_wolfe_variant.use_strong_lazy,
+            linesearch_workspace=workspace,
+            sparsity_control=lazy_tolerance,
+            callback=BDICG_callback,
+            extra_vertex_storage=pre_computed_set,
+        )
+
+    clean_up_pre_computed_set!(lmo, pre_computed_set, x, frank_wolfe_variant)
+
+    return x, primal, dual_gap, status, pre_computed_set
+end
+
+Base.print(io::IO, ::BlendedDecompositionInvariantConditionalGradient) =
+    print(io, "Blended-Decompostion-Invariant-Frank-Wolfe")
+
 """
 	Vanilla-Frank-Wolfe
 

src/interface.jl

@@ -255,7 +255,7 @@ function solve(
         options[:min_fw_iterations],
         time_ref,
         options[:time_limit],
-        use_DICG=typeof(options[:variant]) == DecompositionInvariantConditionalGradient,
+        use_DICG=typeof(options[:variant]) <: DecompositionInvariant,
     )
 
     tree.root.options[:callback] = fw_callback

src/managed_blmo.jl

@@ -106,11 +106,11 @@ function is_inface_feasible(blmo::ManagedBoundedLMO, a, x)
 end
 
 #Provide FrankWolfe.dicg_maximum_step
-function dicg_maximum_step(blmo::ManagedBoundedLMO, x, direction; kwargs...)
+function dicg_maximum_step(blmo::ManagedBoundedLMO, direction, x; kwargs...)
     return bounded_dicg_maximum_step(
         blmo.simple_lmo,
-        x,
         direction,
+        x,
         blmo.lower_bounds,
         blmo.upper_bounds,
         blmo.int_vars,
@@ -218,9 +218,7 @@ end
 # That means does v satisfy all bounds and other linear constraints?
 function is_linear_feasible(blmo::ManagedBoundedLMO, v::AbstractVector)
     for (i, int_var) in enumerate(blmo.int_vars)
-        if !(
-            blmo.lower_bounds[i] ≤ v[int_var] + 1e-6 || !(v[int_var] - 1e-6 ≤ blmo.upper_bounds[i])
-        )
+        if !(blmo.lower_bounds[i] ≤ v[int_var] + 1e-6 && (v[int_var] - 1e-6 ≤ blmo.upper_bounds[i]))
             @debug(
                 "Variable: $(int_var) Vertex entry: $(v[int_var]) Lower bound: $(blmo.lower_bounds[i]) Upper bound: $(blmo.upper_bounds[i]))"
             )

src/node.jl

@@ -133,7 +133,7 @@ function Bonobo.get_branching_nodes_info(tree::Bonobo.BnBTree, node::FrankWolfeN
                 active_set=node.active_set,
                 discarded_vertices=node.discarded_vertices,
                 local_bounds=node.local_bounds,
-                level=node.level + 1,
+                level=(node.level + 1),
                 fw_dual_gap_limit=node.fw_dual_gap_limit,
                 fw_time=node.fw_time,
                 global_tightenings=0,
@@ -158,7 +158,7 @@ function Bonobo.get_branching_nodes_info(tree::Bonobo.BnBTree, node::FrankWolfeN
     end
 
     #different ways to split active set
-    if typeof(tree.root.options[:variant]) != DecompositionInvariantConditionalGradient
+    if !(typeof(tree.root.options[:variant]) <: DecompositionInvariant)
 
         # Keep the same pre_computed_set
         pre_computed_set_left, pre_computed_set_right = node.pre_computed_set, node.pre_computed_set
@@ -181,7 +181,7 @@ function Bonobo.get_branching_nodes_info(tree::Bonobo.BnBTree, node::FrankWolfeN
     discarded_set_left, discarded_set_right =
         split_vertices_set!(node.discarded_vertices, tree, vidx, x, node.local_bounds)
 
-    if typeof(tree.root.options[:variant]) != DecompositionInvariantConditionalGradient
+    if !(typeof(tree.root.options[:variant]) <: DecompositionInvariant)
         # Sanity check
         @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
@@ -227,20 +227,29 @@ 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])
 
-    if typeof(tree.root.options[:variant]) != DecompositionInvariantConditionalGradient
-        # 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)
+    # 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
-        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
-    end
+    active_set_left = build_active_set_by_domain_oracle(
+        x_left,
+        active_set_left,
+        tree,
+        varbounds_left,
+        node;
+        pre_computed_set=pre_computed_set_left,
+        is_decomposition_invariant=typeof(tree.root.options[:variant]) <: DecompositionInvariant,
+    )
+
+    active_set_right = build_active_set_by_domain_oracle(
+        x_right,
+        active_set_right,
+        tree,
+        varbounds_right,
+        node;
+        pre_computed_set=pre_computed_set_right,
+        is_decomposition_invariant=typeof(tree.root.options[:variant]) <: DecompositionInvariant,
+    )
 
     # update the LMO
     node_info_left = (
@@ -343,14 +352,12 @@ function Bonobo.evaluate_node!(tree::Bonobo.BnBTree, node::FrankWolfeNode)
         return NaN, NaN
     end
 
-    if typeof(tree.root.options[:variant]) != DecompositionInvariantConditionalGradient
-        # Check feasibility of the iterate
-        active_set = node.active_set
-        x = FrankWolfe.compute_active_set_iterate!(node.active_set)
-        @assert is_linear_feasible(tree.root.problem.tlmo, x)
-        for (_, v) in node.active_set
-            @assert is_linear_feasible(tree.root.problem.tlmo, v)
-        end
+    # Check feasibility of the iterate
+    active_set = node.active_set
+    x = FrankWolfe.compute_active_set_iterate!(node.active_set)
+    @assert is_linear_feasible(tree.root.problem.tlmo, x)
+    for (_, v) in node.active_set
+        @assert is_linear_feasible(tree.root.problem.tlmo, v)
     end
 
     if tree.root.options[:propagate_bounds] !== nothing
@@ -392,14 +399,6 @@ function Bonobo.evaluate_node!(tree::Bonobo.BnBTree, node::FrankWolfeNode)
     else
         node.pre_computed_set = atoms_set
         node.active_set = FrankWolfe.ActiveSet([(1.0, x)])
-        # update set of computed atoms and active set
-        if isa(x, AbstractVector)
-            node.pre_computed_set = atoms_set
-            node.active_set = FrankWolfe.ActiveSet([(1.0, x)])
-        else
-            @debug "x is not a vector, returning NaN, x: $x"
-            return NaN, NaN
-        end
     end
 
     node.fw_time = Dates.now() - time_ref

src/polytope_blmos.jl

@@ -988,7 +988,7 @@ function is_simple_linear_feasible(sblmo::BirkhoffBLMO, v::AbstractVector)
             return false
         end
         # append by column ? row sum : column sum
-        if !isapprox(sum(v[i:n:n^2]), 1.0, atol=1e-6, rtol=1e-3)
+        if !isapprox(sum(v[i:n:(n^2)]), 1.0, atol=1e-6, rtol=1e-3)
             @debug "Row sum not 1: $(sum(v[i:n:n^2]))"
             return false
         end