diff --git a/Project.toml b/Project.toml index 48761e9db..d3b4bf897 100644 --- a/Project.toml +++ b/Project.toml @@ -22,7 +22,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Bonobo = "0.1.3" DataStructures = "0.18" Distributions = "0.25" -FrankWolfe = "0.2.11" +FrankWolfe = "0.2.24" HiGHS = "1" MathOptInterface = "1" MathOptSetDistances = "0.2" diff --git a/examples/portfolio.jl b/examples/portfolio.jl index 3e09b587b..f5b57741f 100644 --- a/examples/portfolio.jl +++ b/examples/portfolio.jl @@ -53,7 +53,7 @@ const Mi = (Ai + Ai') / 2 return storage end - x, _, result = Boscia.solve(f, grad!, lmo, verbose=true) + x, _, result = Boscia.solve(f, grad!, lmo, verbose=true, time_limit=600) @test dot(ai, x) <= bi + 1e-6 @test f(x) <= f(result[:raw_solution]) + 1e-6 end diff --git a/src/Boscia.jl b/src/Boscia.jl index 08ba19c26..e73ee7b51 100644 --- a/src/Boscia.jl +++ b/src/Boscia.jl @@ -14,6 +14,7 @@ import MathOptSetDistances as MOD include("time_tracking_lmo.jl") include("bounds.jl") +include("frank_wolfe_variants.jl") include("node.jl") include("custom_bonobo.jl") include("callbacks.jl") diff --git a/src/callbacks.jl b/src/callbacks.jl index 629a81f54..03b1d297b 100644 --- a/src/callbacks.jl +++ b/src/callbacks.jl @@ -26,23 +26,25 @@ function build_FW_callback(tree, min_number_lower, check_rounding_value::Bool, f @assert is_linear_feasible(tree.root.problem.lmo, state.v) end # @assert is_linear_feasible(tree.root.problem.lmo, state.x) - push!(fw_iterations, state.t) - if ncalls != state.lmo.ncalls - ncalls = state.lmo.ncalls - (best_v, best_val) = - find_best_solution(tree.root.problem.f, tree.root.problem.lmo.lmo.o, vars) - if best_val < tree.incumbent - tree.root.updated_incumbent[] = true - node = tree.nodes[tree.root.current_node_id[]] - sol = FrankWolfeSolution(best_val, best_v, node, :SCIP) - push!(tree.solutions, sol) - if tree.incumbent_solution === nothing || - sol.objective < tree.incumbent_solution.objective - tree.incumbent_solution = sol + + if state.lmo !== nothing # can happen with using Blended Conditional Gradient + if ncalls != state.lmo.ncalls + ncalls = state.lmo.ncalls + (best_v, best_val) = + find_best_solution(tree.root.problem.f, tree.root.problem.lmo.lmo.o, vars, tree.root.options[:domain_oracle]) + if best_val < tree.incumbent + tree.root.updated_incumbent[] = true + node = tree.nodes[tree.root.current_node_id[]] + sol = FrankWolfeSolution(best_val, best_v, node, :SCIP) + push!(tree.solutions, sol) + if tree.incumbent_solution === nothing || + sol.objective < tree.incumbent_solution.objective + tree.incumbent_solution = sol + end + tree.incumbent = best_val + Bonobo.bound!(tree, node.id) end - tree.incumbent = best_val - Bonobo.bound!(tree, node.id) end end @@ -50,19 +52,21 @@ function build_FW_callback(tree, min_number_lower, check_rounding_value::Bool, f return false end - val = tree.root.problem.f(state.v) - if val < tree.incumbent - tree.root.updated_incumbent[] = true - #TODO: update solution without adding node - node = tree.nodes[tree.root.current_node_id[]] - sol = FrankWolfeSolution(val, copy(state.v), node, :vertex) - push!(tree.solutions, sol) - if tree.incumbent_solution === nothing || - sol.objective < tree.incumbent_solution.objective - tree.incumbent_solution = sol + if tree.root.options[:domain_oracle](state.v) + val = tree.root.problem.f(state.v) + if val < tree.incumbent + tree.root.updated_incumbent[] = true + #TODO: update solution without adding node + node = tree.nodes[tree.root.current_node_id[]] + sol = FrankWolfeSolution(val, copy(state.v), node, :vertex) + 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 - tree.incumbent = val - Bonobo.bound!(tree, node.id) end node = tree.nodes[tree.root.current_node_id[]] diff --git a/src/custom_bonobo.jl b/src/custom_bonobo.jl index 837378bd8..0fcf91775 100644 --- a/src/custom_bonobo.jl +++ b/src/custom_bonobo.jl @@ -106,5 +106,12 @@ function Bonobo.get_solution( tree::Bonobo.BnBTree{N,R,V,S}; result=1, ) where {N,R,V,S<:FrankWolfeSolution{N,V}} + if isempty(tree.solutions) + @warn "There is no solution in the tree. This behaviour can happen if you have supplied + \na custom domain oracle. In that case, try to increase the time limit. If you have not specified a + \ndomain oracle, please report!" + @assert tree.root.problem.solving_stage == TIME_LIMIT_REACHED + return nothing + end return tree.solutions[result].solution end diff --git a/src/frank_wolfe_variants.jl b/src/frank_wolfe_variants.jl new file mode 100644 index 000000000..0c18a6318 --- /dev/null +++ b/src/frank_wolfe_variants.jl @@ -0,0 +1,217 @@ + +""" +Frank-Wolfe variant used to compute the problems at node level. +A `FrankWolfeVariant` must implement +``` +solve_frank_wolfe(fw::FrankWolfeVariant, f, grad!, lmo, active_set, line_search, epsilon, max_iteration, + added_dropped_vertices, use_extra_vertex_storage, callback, lazy, timeout, verbose, workspace)) +``` +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 + +# default printing for FrankWolfeVariant is just showing the type +Base.print(io::IO, fw::FrankWolfeVariant) = print(io, split(string(typeof(fw)), ".")[end]) + +""" + solve_frank_wolfe(fw::FrankWolfeVariant, f, grad!, lmo, active_set, line_search, epsilon, max_iteration, + added_dropped_vertices, use_extra_vertex_storage, callback, lazy, timeout, verbose, workspace) + +Returns the optimal solution x to the node problem, its primal and dual gap and the active set. +""" +function solve_frank_wolfe end + +build_frank_wolfe_workspace(::FrankWolfeVariant, x) = nothing + + +""" + Away-Frank-Wolfe + +In every iteration, it computes the worst performing vertex, called away vertex, in the active set with regard to the gradient. +If enough local progress can be made, weight is shifted from the away vertex to all other vertices. + +In case lazification is activated, the FW vertex is only computed if not enough local progress can be guaranteed. +""" +struct AwayFrankWolfe <: FrankWolfeVariant end + +function solve_frank_wolfe( + frank_wolfe_variant::AwayFrankWolfe, + f, + grad!, + lmo, + active_set; + line_search::FrankWolfe.LineSearchMethod=FrankWolfe.Adaptive(), + epsilon=1e-7, + max_iteration=10000, + add_dropped_vertices=false, + use_extra_vertex_storage=false, + extra_vertex_storage=nothing, + callback=nothing, + lazy=false, + timeout=Inf, + verbose=false, + workspace=nothing +) + x, _, primal, dual_gap, _, active_set = FrankWolfe.away_frank_wolfe( + f, + grad!, + lmo, + active_set, + epsilon=epsilon, + max_iteration=max_iteration, + line_search=line_search, + callback=callback, + lazy=lazy, + timeout=timeout, + add_dropped_vertices=add_dropped_vertices, + use_extra_vertex_storage=use_extra_vertex_storage, + extra_vertex_storage=extra_vertex_storage, + verbose=verbose, + ) + + return x, primal, dual_gap, active_set +end + +Base.print(io::IO, ::AwayFrankWolfe) = print(io, "Away-Frank-Wolfe") + + +""" + Blended Conditional Gradient +""" +struct Blended <: FrankWolfeVariant end + +function solve_frank_wolfe( + frank_wolfe_variant::Blended, + f, + grad!, + lmo, + active_set; + line_search::FrankWolfe.LineSearchMethod=FrankWolfe.Adaptive(), + epsilon=1e-7, + max_iteration=10000, + add_dropped_vertices=false, + use_extra_vertex_storage=false, + extra_vertex_storage=nothing, + callback=nothing, + lazy=false, + timeout=Inf, + verbose=false, + workspace=nothing +) + x,_, primal, dual_gap,_, active_set = blended_conditional_gradient( + f, + grad!, + lmo, + active_set, + line_search=line_search, + epsilon=epsilon, + max_iteration=max_iteration, + add_dropped_vertices=add_dropped_vertices, + use_extra_vertex_storage=use_extra_vertex_storage, + extra_vertex_storage=extra_vertex_storage, + callback=callback, + timeout=timeout, + verbose=verbose, + ) + + return x, primal, dual_gap, active_set +end + +Base.print(io::IO, ::Blended) = print(io, "Blended Conditional Gradient") + +""" + Blended Pairwise Conditional Gradient +""" +struct BPCG <: FrankWolfeVariant end + +function solve_frank_wolfe( + frank_wolfe_variant::BPCG, + f, + grad!, + lmo, + active_set; + line_search::FrankWolfe.LineSearchMethod=FrankWolfe.Adaptive(), + epsilon=1e-7, + max_iteration=10000, + add_dropped_vertices=false, + use_extra_vertex_storage=false, + extra_vertex_storage=nothing, + callback=nothing, + lazy=false, + timeout=Inf, + verbose=false, + workspace=nothing +) + x, _, primal, dual_gap, _, active_set = FrankWolfe.blended_pairwise_conditional_gradient( + f, + grad!, + lmo, + active_set, + line_search=line_search, + epsilon=epsilon, + max_iteration=max_iteration, + add_dropped_vertices=add_dropped_vertices, + use_extra_vertex_storage=use_extra_vertex_storage, + extra_vertex_storage=extra_vertex_storage, + callback=callback, + lazy=lazy, + timeout=timeout, + verbose=verbose + ) + + return x, primal, dual_gap, active_set +end + +Base.print(io::IO, ::BPCG) = print(io, "Blended Pairwise Conditional Gradient") + +""" + Vanilla-Frank-Wolfe + +The standard variant of Frank-Wolfe. In each iteration, the vertex v minimizing ∇f * (x-v) is computed. + +Lazification cannot be used in this setting. +""" +struct VanillaFrankWolfe <: FrankWolfeVariant end + +function solve_frank_wolfe( + frank_wolfe_variant::VanillaFrankWolfe, + f, + grad!, + lmo, + active_set; + line_search::FrankWolfe.LineSearchMethod=FrankWolfe.Adaptive(), + epsilon=1e-7, + max_iteration=10000, + add_dropped_vertices=false, + use_extra_vertex_storage=false, + extra_vertex_storage=nothing, + callback=nothing, + lazy=false, + timeout=Inf, + verbose=false, + workspace=nothing +) + # If the flag away_steps is set to false, away_frank_wolfe performs Vanilla. + # Observe that the lazy flag is only observed if away_steps is set to true, so it can neglected. + x, _, primal, dual_gap, _, active_set = FrankWolfe.away_frank_wolfe( + f, + grad!, + lmo, + active_set, + away_steps=false, + epsilon=epsilon, + max_iteration=max_iteration, + line_search=line_search, + callback=callback, + timeout=timeout, + add_dropped_vertices=add_dropped_vertices, + use_extra_vertex_storage=use_extra_vertex_storage, + extra_vertex_storage=extra_vertex_storage, + verbose=verbose, + ) + + return x, primal, dual_gap, active_set +end + +Base.print(io::IO, ::VanillaFrankWolfe) = print(io, "Vanilla-Frank-Wolfe") diff --git a/src/heuristics.jl b/src/heuristics.jl index f9e412b34..ebbcd47a4 100644 --- a/src/heuristics.jl +++ b/src/heuristics.jl @@ -3,34 +3,38 @@ Finds the best solution in the SCIP solution storage, based on the objective function `f`. Returns the solution vector and the corresponding best value. """ -function find_best_solution(f::Function, o::SCIP.Optimizer, vars::Vector{MOI.VariableIndex}) +function find_best_solution(f::Function, o::SCIP.Optimizer, vars::Vector{MOI.VariableIndex}, domain_oracle) sols_vec = unsafe_wrap(Vector{Ptr{Cvoid}}, SCIP.LibSCIP.SCIPgetSols(o), SCIP.LibSCIP.SCIPgetNSols(o)) best_val = Inf best_v = nothing for sol in sols_vec v = SCIP.sol_values(o, vars, sol) - val = f(v) - if val < best_val - best_val = val - best_v = v + if domain_oracle(v) + val = f(v) + if val < best_val + best_val = val + best_v = v + end end end - @assert isfinite(best_val) + #@assert isfinite(best_val) -> not necessarily the case if the domain oracle is not the default. return (best_v, best_val) end -function find_best_solution(f::Function, o::MOI.AbstractOptimizer, vars::Vector{MOI.VariableIndex}) +function find_best_solution(f::Function, o::MOI.AbstractOptimizer, vars::Vector{MOI.VariableIndex}, domain_oracle) nsols = MOI.get(o, MOI.ResultCount()) @assert nsols > 0 best_val = Inf best_v = nothing for sol_idx in 1:nsols xv = [MOI.get(o, MOI.VariablePrimal(sol_idx), xi) for xi in vars] - val = f(xv) - if val < best_val - best_val = val - best_v = xv + if domain_oracle(xv) + val = f(xv) + if val < best_val + best_val = val + best_v = xv + end end end return (best_v, best_val) diff --git a/src/interface.jl b/src/interface.jl index a30d1a39e..07bab846e 100644 --- a/src/interface.jl +++ b/src/interface.jl @@ -1,34 +1,44 @@ """ solve -f - objective function oracle. -g - oracle for the gradient of the objective. -lmo - a MIP solver instance (SCIP) encoding the feasible region. -traverse_strategy - encodes how to choose the next node for evaluation. +f - objective function oracle. +g - oracle for the gradient of the objective. +lmo - a MIP solver instance (SCIP) encoding the feasible region. +traverse_strategy - encodes how to choose the next node for evaluation. By default the node with the best lower bound is picked. -branching_strategy - by default we branch on the entry which is the farthest +branching_strategy - by default we branch on the entry which is the farthest away from being an integer. -fw_epsilon - the tolerance for FrankWolfe in the root node. -verbose - if true, a log and solution statistics are printed. -dual_gap - if this absolute dual gap is reached, the algorithm stops. -rel_dual_gap - if this relative dual gap is reached, the algorithm stops. -time_limit - algorithm will stop if the time limit is reached. Depending on the problem +variant - variant of FrankWolfe to be used to solve the node problem. + Options: FW -- Vanilla FrankWolfe + AFW -- Away FrankWolfe + BPCG -- Blended Pairwise Conditional Gradient +line_search - specifies the Line Search method used in the FrankWolfe variant. + Default is the Adaptive Line Search. For other types, check the FrankWolfe.jl package. +fw_epsilon - the tolerance for FrankWolfe in the root node. +verbose - if true, a log and solution statistics are printed. +dual_gap - if this absolute dual gap is reached, the algorithm stops. +rel_dual_gap - if this relative dual gap is reached, the algorithm stops. +time_limit - algorithm will stop if the time limit is reached. Depending on the problem it is possible that no feasible solution has been found yet. -print_iter - encodes after how many proccessed nodes the current node and solution status +print_iter - encodes after how many proccessed nodes the current node and solution status is printed. Will always print if a new integral solution has been found. -dual_gap_decay_factor - the FrankWolfe tolerance at a given level i in the tree is given by +dual_gap_decay_factor - the FrankWolfe tolerance at a given level i in the tree is given by fw_epsilon * dual_gap_decay_factor^i until we reach the min_node_fw_epsilon. -max_fw_iter - maximum number of iterations in a FrankWolfe run. -min_number_lower - If not Inf, evaluation of a node is stopped if at least min_number_lower nodes have a better +max_fw_iter - maximum number of iterations in a FrankWolfe run. +min_number_lower - If not Inf, evaluation of a node is stopped if at least min_number_lower nodes have a better lower bound. -min_node_fw_epsilon - smallest fw epsilon possible, see dual_gap_decay_factor. -min_fw_iterations - the minimum number of FrankWolfe iterations in the node evaluation. -max_iteration_post - maximum number of iterations in a FrankWolfe run during postsolve -dual_tightening - whether to use dual tightening techniques (make sure your function is convex!) +min_node_fw_epsilon - smallest fw epsilon possible, see dual_gap_decay_factor. +min_fw_iterations - the minimum number of FrankWolfe iterations in the node evaluation. +max_iteration_post - maximum number of iterations in a FrankWolfe run during postsolve +dual_tightening - whether to use dual tightening techniques (make sure your function is convex!) global_dual_tightening - dual tightening maintained globally valid (when new solutions are found) -bnb_callback - an optional callback called at every node of the tree, for example for heuristics -strong_convexity - strong convexity of the function, used for tighter dual bound at every node -start_solution - initial solution to start with an incumbent +bnb_callback - an optional callback called at every node of the tree, for example for heuristics +strong_convexity - strong convexity of the function, used for tighter dual bound at every node +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 + on the Line Search method, you might have to provide the domain oracle to it, too. +start_solution - initial solution to start with an incumbent +fw_verbose - if true, FrankWolfe logs are printed """ function solve( f, @@ -36,6 +46,9 @@ function solve( lmo; traverse_strategy=Bonobo.BestFirstSearch(), branching_strategy=Bonobo.MOST_INFEASIBLE(), + variant::FrankWolfeVariant=BPCG(), + line_search::FrankWolfe.LineSearchMethod=FrankWolfe.Adaptive(), + active_set::Union{Nothing, FrankWolfe.ActiveSet} = nothing, fw_epsilon=1e-2, verbose=false, dual_gap=1e-6, @@ -53,7 +66,9 @@ function solve( global_dual_tightening=true, bnb_callback=nothing, strong_convexity=0.0, + domain_oracle= x->true, start_solution=nothing, + fw_verbose = false, kwargs..., ) if verbose @@ -61,6 +76,8 @@ function solve( println("Parameter settings.") println("\t Tree traversal strategy: ", _value_to_print(traverse_strategy)) println("\t Branching strategy: ", _value_to_print(branching_strategy)) + println("\t FrankWolfe variant: $(variant)") + println("\t Line Search Method: $(line_search)") @printf("\t Absolute dual gap tolerance: %e\n", dual_gap) @printf("\t Relative dual gap tolerance: %e\n", rel_dual_gap) @printf("\t Frank-Wolfe subproblem tolerance: %e\n", fw_epsilon) @@ -136,10 +153,20 @@ function solve( global_bounds[idx, :lessthan] = MOI.LessThan(1.0) end - direction = collect(1.0:n) - v = compute_extreme_point(lmo, direction) - v[integer_variables] = round.(v[integer_variables]) - active_set = FrankWolfe.ActiveSet([(1.0, v)]) + v = [] + if active_set === nothing + direction = collect(1.0:n) + v = compute_extreme_point(lmo, direction) + v[integer_variables] = round.(v[integer_variables]) + active_set = FrankWolfe.ActiveSet([(1.0, v)]) + vertex_storage = FrankWolfe.DeletedVertexStorage(typeof(v)[], 1) + else + @assert FrankWolfe.active_set_validate(active_set) + for a in active_set.atoms + @assert is_linear_feasible(lmo.o, a) + end + v = active_set.atoms[1] + end vertex_storage = FrankWolfe.DeletedVertexStorage(typeof(v)[], 1) m = Boscia.SimpleOptimizationProblem(f, grad!, n, integer_variables, time_lmo, global_bounds) @@ -183,6 +210,11 @@ function solve( :dual_tightening => dual_tightening, :global_dual_tightening => global_dual_tightening, :strong_convexity => strong_convexity, + :variant => variant, + :lineSearch => line_search, + :domain_oracle => domain_oracle, + :usePostsolve => use_postsolve, + :fwVerbose => fw_verbose, ), ), branch_strategy=branching_strategy, @@ -269,22 +301,24 @@ function solve( Bonobo.optimize!(tree; callback=bnb_callback) - x = postsolve(tree, result, time_ref, verbose, use_postsolve, max_iteration_post) + x = postsolve(tree, result, time_ref, verbose, max_iteration_post) # Check solution and polish x_polished = x - if !is_linear_feasible(tree.root.problem.lmo, x) - error("Reported solution not linear feasbile!") - end - if !is_integer_feasible(tree.root.problem.integer_variables, x, atol=1e-16, rtol=1e-16) - @info "Polish solution" - for i in tree.root.problem.integer_variables - x_polished[i] = round(x_polished[i]) + if x !== nothing + if !is_linear_feasible(tree.root.problem.lmo, x) + error("Reported solution not linear feasbile!") end - if !is_linear_feasible(tree.root.problem.lmo, x_polished) - @warn "Polished solution not linear feasible" - else - x = x_polished + if !is_integer_feasible(tree.root.problem.integer_variables, x, atol=1e-16, rtol=1e-16) && x !== nothing + @info "Polish solution" + for i in tree.root.problem.integer_variables + x_polished[i] = round(x_polished[i]) + end + if !is_linear_feasible(tree.root.problem.lmo, x_polished) + @warn "Polished solution not linear feasible" + else + x = x_polished + end end end println() # cleaner output @@ -492,7 +526,13 @@ function build_bnb_callback( if Bonobo.terminated(tree) Bonobo.sort_solutions!(tree.solutions, tree.sense) x = Bonobo.get_solution(tree) - primal_value = tree.root.problem.f(x) + # x can be nothing if the user supplied a custom domain oracle and the time limit is reached + if x === nothing + @assert tree.root.problem.solving_stage == TIME_LIMIT_REACHED + end + primal_value = x !== nothing ? tree.root.problem.f(x) : Inf + # deactivate postsolve if there is no solution + tree.root.options[:usePostsolve] = x === nothing ? false : tree.root.options[:usePostsolve] # TODO: here we need to calculate the actual state @@ -533,9 +573,9 @@ Runs the post solve both for a cleaner solutiona and to optimize for the continuous variables if present. Prints solution statistics if verbose is true. """ -function postsolve(tree, result, time_ref, verbose, use_postsolve, max_iteration_post) +function postsolve(tree, result, time_ref, verbose, max_iteration_post) x = Bonobo.get_solution(tree) - primal = tree.incumbent_solution.objective + primal = x !== nothing ? tree.incumbent_solution.objective : Inf status_string = "FIX ME" # should report "feasible", "optimal", "infeasible", "gap tolerance met" if isempty(tree.nodes) @@ -549,7 +589,7 @@ function postsolve(tree, result, time_ref, verbose, use_postsolve, max_iteration end only_integer_vars = tree.root.problem.nvars == length(tree.root.problem.integer_variables) - if use_postsolve && !only_integer_vars + if tree.root.options[:usePostsolve] && !only_integer_vars # Build solution lmo fix_bounds = IntegerBounds() for i in tree.root.problem.integer_variables diff --git a/src/node.jl b/src/node.jl index f066bda3f..817200a92 100644 --- a/src/node.jl +++ b/src/node.jl @@ -72,11 +72,11 @@ function Bonobo.get_branching_nodes_info(tree::Bonobo.BnBTree, node::FrankWolfeN end new_bound_left = lower_bound_base + μ/2 * (x[vidx] - floor(x[vidx]))^2 new_bound_right = lower_bound_base + μ/2 * (ceil(x[vidx]) - x[vidx])^2 - if new_bound_left >= tree.incumbent + if new_bound_left > tree.incumbent @debug "prune left, from $(node.lb) -> $new_bound_left, ub $(tree.incumbent), lb $(node.lb)" prune_left = true end - if new_bound_right >= tree.incumbent + if new_bound_right > tree.incumbent @debug "prune right, from $(node.lb) -> $new_bound_right, ub $(tree.incumbent), lb $(node.lb)" prune_right = true end @@ -157,11 +157,21 @@ function Bonobo.get_branching_nodes_info(tree::Bonobo.BnBTree, node::FrankWolfeN local_potential_tightenings=0, dual_gap=NaN, ) - nodes = if !prune_left && !prune_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) + domain_oracle = tree.root.options[:domain_oracle] + + nodes = if !prune_left && !prune_right #&& domain_oracle(x_left) && domain_oracle(x_right) [node_info_left, node_info_right] - elseif prune_left + elseif prune_left [node_info_right] - else + elseif prune_right + [node_info_left] + elseif domain_oracle(x_right) + [node_info_right] + elseif domain_oracle(x_left) [node_info_left] end return nodes @@ -262,21 +272,63 @@ function Bonobo.evaluate_node!(tree::Bonobo.BnBTree, node::FrankWolfeNode) end end - # call blended_pairwise_conditional_gradient - x, _, primal, dual_gap, _, active_set = FrankWolfe.blended_pairwise_conditional_gradient( + # x = zeros(tree.root.problem.nvars) + #dual_gap = primal = 0 + #active_set = FrankWolfe.ActiveSet([(1.0, x)]) + x=FrankWolfe.compute_active_set_iterate!(node.active_set) + domain_oracle = tree.root.options[:domain_oracle] + + x, primal, dual_gap, active_set = solve_frank_wolfe( + tree.root.options[:variant], tree.root.problem.f, tree.root.problem.g, tree.root.problem.lmo, - node.active_set, + node.active_set; epsilon=node.fw_dual_gap_limit, max_iteration=tree.root.options[:max_fw_iter], - line_search=FrankWolfe.Adaptive(verbose=false), + line_search=tree.root.options[:lineSearch], + lazy=true, + timeout=tree.root.options[:time_limit], add_dropped_vertices=true, use_extra_vertex_storage=true, extra_vertex_storage=node.discarded_vertices, callback=tree.root.options[:callback], - lazy=true, - ) + verbose=tree.root.options[:fwVerbose], + )#= + if tree.root.options[:variant] == BPCG + # call blended_pairwise_conditional_gradient + x, _, primal, dual_gap, _, active_set = FrankWolfe.blended_pairwise_conditional_gradient( + tree.root.problem.f, + tree.root.problem.g, + tree.root.problem.lmo, + node.active_set, + epsilon=node.fw_dual_gap_limit, + max_iteration=tree.root.options[:max_fw_iter], + line_search=tree.root.options[:lineSearch], + add_dropped_vertices=true, + use_extra_vertex_storage=true, + extra_vertex_storage=node.discarded_vertices, + callback=tree.root.options[:callback], + lazy=true, + timeout=tree.root.options[:time_limit], + verbose=tree.root.options[:fwVerbose], + ) + elseif tree.root.options[:variant] == AFW + # call away_frank_wolfe + x, _, primal, dual_gap, _, active_set = FrankWolfe.away_frank_wolfe( + tree.root.problem.f, + tree.root.problem.g, + tree.root.problem.lmo, + node.active_set, + epsilon=node.fw_dual_gap_limit, + max_iteration=tree.root.options[:max_fw_iter], + line_search=tree.root.options[:lineSearch], + callback=tree.root.options[:callback], + lazy=true, + timeout=tree.root.options[:time_limit], + verbose=tree.root.options[:fwVerbose], + ) + end=# node.fw_time = Dates.now() - time_ref node.dual_gap = dual_gap diff --git a/test/interface_test.jl b/test/interface_test.jl index 867fa2ad0..d3b5039f0 100644 --- a/test/interface_test.jl +++ b/test/interface_test.jl @@ -2,6 +2,7 @@ using LinearAlgebra using Distributions import Random using SCIP +using HiGHS import MathOptInterface const MOI = MathOptInterface import Boscia @@ -170,3 +171,107 @@ Ns = 0.1 @test sum(x2[p+1:2p]) <= k @test f(x2) == f(x) end + +n = 20 +diffi = Random.rand(Bool, n) * 0.6 .+ 0.3 + +@testset "Different FW variants" begin + + function build_model() + o = SCIP.Optimizer() + MOI.set(o, MOI.Silent(), true) + MOI.empty!(o) + x = MOI.add_variables(o, n) + for xi in x + MOI.add_constraint(o, xi, MOI.GreaterThan(0.0)) + MOI.add_constraint(o, xi, MOI.LessThan(1.0)) + MOI.add_constraint(o, xi, MOI.ZeroOne()) # or MOI.Integer() + end + lmo = FrankWolfe.MathOptLMO(o) + return lmo + end + + 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 + + direction = rand(n) + + lmo = build_model() + v = FrankWolfe.compute_extreme_point(lmo, direction) + x_afw, _, result_afw = Boscia.solve(f, grad!, lmo, verbose=false, variant=Boscia.AwayFrankWolfe()) + + lmo = build_model() + v = FrankWolfe.compute_extreme_point(lmo, direction) + x_blended, _, result_blended = Boscia.solve(f, grad!, lmo, verbose=false, variant=Boscia.Blended()) + + lmo = build_model() + v = FrankWolfe.compute_extreme_point(lmo, direction) + x_bpcg, _, result_bpcg = Boscia.solve(f, grad!, lmo, verbose=false, variant=Boscia.BPCG()) + + lmo = build_model() + v = FrankWolfe.compute_extreme_point(lmo, direction) + x_vfw, _, result_vfw = Boscia.solve(f, grad!, lmo, verbose=false, variant=Boscia.VanillaFrankWolfe()) + + @test isapprox(f(x_afw), f(result_afw[:raw_solution]), atol=1e-6, rtol=1e-3) + @test isapprox(f(x_blended), f(result_blended[:raw_solution]), atol=1e-6, rtol=1e-3) + @test isapprox(f(x_bpcg), f(result_bpcg[:raw_solution]), atol=1e-6, rtol=1e-3) + @test isapprox(f(x_vfw), f(result_vfw[:raw_solution]), atol=1e-6, rtol=1e-3) + + @test sum(isapprox.(x_afw, x_blended, atol=1e-6, rtol=1e-3)) == n + @test sum(isapprox.(x_blended, x_bpcg, atol=1e-6, rtol=1e-3)) == n + @test sum(isapprox.(x_bpcg, x_vfw, atol=1e-6, rtol=1e-3)) == n + @test sum(isapprox.(x_vfw, x_afw, atol=1e-6, rtol=1e-3)) == n +end + +@testset "Different line search types" begin + + function build_model() + o = SCIP.Optimizer() + MOI.set(o, MOI.Silent(), true) + MOI.empty!(o) + x = MOI.add_variables(o, n) + for xi in x + MOI.add_constraint(o, xi, MOI.GreaterThan(0.0)) + MOI.add_constraint(o, xi, MOI.LessThan(1.0)) + MOI.add_constraint(o, xi, MOI.ZeroOne()) # or MOI.Integer() + end + lmo = FrankWolfe.MathOptLMO(o) + return lmo + end + + 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 + + direction = rand(n) + + lmo = build_model() + v = FrankWolfe.compute_extreme_point(lmo, direction) + line_search = FrankWolfe.Adaptive() + x_adaptive, _, result_adaptive = Boscia.solve(f, grad!, lmo, verbose=false, line_search=line_search) + + lmo = build_model() + v = FrankWolfe.compute_extreme_point(lmo, direction) + line_search = FrankWolfe.MonotonicStepSize() + x_monotonic, _, result_monotonic = Boscia.solve(f, grad!, lmo, verbose=false, line_search=line_search, time_limit=120) + + lmo = build_model() + v = FrankWolfe.compute_extreme_point(lmo, direction) + line_search = FrankWolfe.Agnostic() + x_agnostic, _, result_agnostic = Boscia.solve(f, grad!, lmo, verbose=false, line_search=line_search, time_limit=120) + + @test isapprox(f(x_adaptive), f(result_adaptive[:raw_solution]), atol=1e-6, rtol=1e-3) + @test isapprox(f(x_monotonic), f(result_monotonic[:raw_solution]), atol=1e-6, rtol=1e-3) + @test isapprox(f(x_agnostic), f(result_agnostic[:raw_solution]), atol=1e-6, rtol=1e-3) + + @test sum(isapprox.(x_adaptive, x_monotonic, atol=1e-6, rtol=1e-3)) == n + @test sum(isapprox.(x_agnostic, x_monotonic, atol=1e-6, rtol=1e-3)) == n + @test sum(isapprox.(x_adaptive, x_agnostic, atol=1e-6, rtol=1e-3)) == n +end \ No newline at end of file diff --git a/test/mean_risk.jl b/test/mean_risk.jl index 178e653c9..c51eb53b6 100644 --- a/test/mean_risk.jl +++ b/test/mean_risk.jl @@ -61,6 +61,6 @@ const M1 = (A1 + A1') / 2 x,_,result = Boscia.solve(f, grad!, lmo, verbose= true) - @test dot(a,x) <= b + 1e-6 + @test dot(a,x) <= b + 1e-4 @test f(x) <= f(result[:raw_solution]) + 1e-6 end diff --git a/test/runtests.jl b/test/runtests.jl index d3e0966f7..df0d92e79 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -131,6 +131,7 @@ end @test isapprox(x, round.(diff1), atol=1e-5, rtol=1e-5) end + for file in readdir(joinpath(@__DIR__, "../examples/"), join=true) if endswith(file, "jl") include(file)