From 367e73ff122a3cda507f5b1ead6013ed71b7166c Mon Sep 17 00:00:00 2001 From: Hendrych Date: Fri, 30 Aug 2024 12:59:37 +0200 Subject: [PATCH 01/45] Move tightening functions into a separate file. --- src/Boscia.jl | 1 + src/node.jl | 249 --------------------------------------------- src/tightenings.jl | 247 ++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 248 insertions(+), 249 deletions(-) create mode 100644 src/tightenings.jl diff --git a/src/Boscia.jl b/src/Boscia.jl index 5e92fc6ca..6853fa541 100644 --- a/src/Boscia.jl +++ b/src/Boscia.jl @@ -17,6 +17,7 @@ include("blmo_interface.jl") include("time_tracking_lmo.jl") include("frank_wolfe_variants.jl") include("build_lmo.jl") +include("tightenings.jl") include("node.jl") include("custom_bonobo.jl") include("callbacks.jl") diff --git a/src/node.jl b/src/node.jl index c8af87ac3..66fe189e1 100644 --- a/src/node.jl +++ b/src/node.jl @@ -177,47 +177,6 @@ function Bonobo.get_branching_nodes_info(tree::Bonobo.BnBTree, node::FrankWolfeN return nodes end -""" -Use strong convexity to potentially remove one of the children nodes -""" -function prune_children(tree, node, lower_bound_base, x, vidx) - prune_left = false - prune_right = false - - μ = tree.root.options[:strong_convexity] - if μ > 0 - @debug "Using strong convexity $μ" - for j in tree.root.problem.integer_variables - if vidx == j - continue - end - lower_bound_base += μ / 2 * min((x[j] - floor(x[j]))^2, (ceil(x[j]) - x[j])^2) - 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 - @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 - @debug "prune right, from $(node.lb) -> $new_bound_right, ub $(tree.incumbent), lb $(node.lb)" - prune_right = true - end - @assert !( - (new_bound_left > tree.incumbent + tree.root.options[:dual_gap]) && - (new_bound_right > tree.incumbent + tree.root.options[:dual_gap]) - ) "both sides should not be pruned" - end - - # If both nodes are pruned, when one of them has to be equal to the incumbent. - # Thus, we have proof of optimality by strong convexity. - if prune_left && prune_right - tree.lb = min(new_bound_left, new_bound_right) - end - - return prune_left, prune_right -end - """ Computes the relaxation at that node """ @@ -329,214 +288,6 @@ function Bonobo.evaluate_node!(tree::Bonobo.BnBTree, node::FrankWolfeNode) return lower_bound, NaN end -""" -Tightening of the bounds at node level. Children node inherit the updated bounds. -""" -function dual_tightening(tree, node, x, dual_gap) - if tree.root.options[:dual_tightening] && isfinite(tree.incumbent) - grad = similar(x) - tree.root.problem.g(grad, x) - num_tightenings = 0 - num_potential_tightenings = 0 - μ = tree.root.options[:strong_convexity] - for j in tree.root.problem.integer_variables - lb_global = get(tree.root.problem.integer_variable_bounds, (j, :greaterthan), -Inf) - ub_global = get(tree.root.problem.integer_variable_bounds, (j, :lessthan), Inf) - lb = get(node.local_bounds.lower_bounds, j, lb_global) - ub = get(node.local_bounds.upper_bounds, j, ub_global) - @assert lb >= lb_global - @assert ub <= ub_global - if lb ≈ ub - # variable already fixed - continue - end - gj = grad[j] - safety_tolerance = 2.0 - rhs = - tree.incumbent - tree.root.problem.f(x) + - safety_tolerance * dual_gap + - sqrt(eps(tree.incumbent)) - if ≈(x[j], lb, atol=tree.options.atol, rtol=tree.options.rtol) - if !isapprox(gj, 0, atol=1e-5) - num_potential_tightenings += 1 - end - if gj > 0 - Mlb = 0 - bound_tightened = true - @debug "starting tightening ub $(rhs)" - while 0.99 * (Mlb * gj + μ / 2 * Mlb^2) <= rhs - Mlb += 1 - if lb + Mlb - 1 == ub - bound_tightened = false - break - end - end - if bound_tightened - new_bound = lb + Mlb - 1 - @debug "found UB tightening $ub -> $new_bound" - node.local_bounds[j, :lessthan] = new_bound - num_tightenings += 1 - if haskey(tree.root.problem.integer_variable_bounds, (j, :lessthan)) - @assert node.local_bounds[j, :lessthan] <= - tree.root.problem.integer_variable_bounds[j, :lessthan] - end - end - end - elseif ≈(x[j], ub, atol=tree.options.atol, rtol=tree.options.rtol) - if !isapprox(gj, 0, atol=1e-5) - num_potential_tightenings += 1 - end - if gj < 0 - Mub = 0 - bound_tightened = true - @debug "starting tightening lb $(rhs)" - while -0.99 * (Mub * gj + μ / 2 * Mub^2) <= rhs - Mub += 1 - if ub - Mub + 1 == lb - bound_tightened = false - break - end - end - if bound_tightened - new_bound = ub - Mub + 1 - @debug "found LB tightening $lb -> $new_bound" - node.local_bounds[j, :greaterthan] = new_bound - num_tightenings += 1 - if haskey(tree.root.problem.integer_variable_bounds, (j, :greaterthan)) - @assert node.local_bounds[j, :greaterthan] >= - tree.root.problem.integer_variable_bounds[j, :greaterthan] - end - end - end - end - end - @debug "# tightenings $num_tightenings" - node.local_tightenings = num_tightenings - node.local_potential_tightenings = num_potential_tightenings - end -end - -""" -Save the gradient of the root solution (i.e. the relaxed solution) and the -corresponding lower and upper bounds. -""" -function store_data_global_tightening(tree, node, x, dual_gap) - if tree.root.options[:global_dual_tightening] && node.std.id == 1 - @debug "storing root node info for tightening" - grad = similar(x) - tree.root.problem.g(grad, x) - safety_tolerance = 2.0 - tree.root.global_tightening_rhs[] = -tree.root.problem.f(x) + safety_tolerance * dual_gap - for j in tree.root.problem.integer_variables - if haskey(tree.root.problem.integer_variable_bounds.upper_bounds, j) - ub = tree.root.problem.integer_variable_bounds[j, :lessthan] - if ≈(x[j], ub, atol=tree.options.atol, rtol=tree.options.rtol) && grad[j] < 0 - tree.root.global_tightening_root_info.upper_bounds[j] = (grad[j], ub) - end - end - if haskey(tree.root.problem.integer_variable_bounds.lower_bounds, j) - lb = tree.root.problem.integer_variable_bounds[j, :greaterthan] - if ≈(x[j], lb, atol=tree.options.atol, rtol=tree.options.rtol) && grad[j] > 0 - tree.root.global_tightening_root_info.lower_bounds[j] = (grad[j], lb) - end - end - end - end -end - -""" -Use the gradient of the root node to tighten the global bounds. -""" -function global_tightening(tree, node) - # new incumbent: check global fixings - if tree.root.options[:global_dual_tightening] && tree.root.updated_incumbent[] - num_tightenings = 0 - rhs = tree.incumbent + tree.root.global_tightening_rhs[] - @assert isfinite(rhs) - for (j, (gj, lb)) in tree.root.global_tightening_root_info.lower_bounds - ub = get(tree.root.problem.integer_variable_bounds.upper_bounds, j, Inf) - ub_new = get(tree.root.global_tightening_root_info.upper_bounds, j, Inf) - ub = min(ub, ub_new) - Mlb = 0 - bound_tightened = true - lb = lb - while Mlb * gj <= rhs - Mlb += 1 - if lb + Mlb - 1 == ub - bound_tightened = false - break - end - end - if bound_tightened - new_bound = lb + Mlb - 1 - @debug "found global UB tightening $ub -> $new_bound" - if haskey(tree.root.global_tightenings.upper_bounds, j) - if tree.root.global_tightenings.upper_bounds[j] != new_bound - num_tightenings += 1 - end - else - num_tightenings += 1 - end - tree.root.global_tightenings.upper_bounds[j] = new_bound - end - end - for (j, (gj, ub)) in tree.root.global_tightening_root_info.upper_bounds - lb = get(tree.root.problem.integer_variable_bounds.lower_bounds, j, -Inf) - lb_new = get(tree.root.global_tightening_root_info.lower_bounds, j, -Inf) - lb = max(lb, lb_new) - Mub = 0 - bound_tightened = true - ub = ub - while -Mub * gj <= rhs - Mub += 1 - if ub - Mub + 1 == lb - bound_tightened = false - break - end - end - if bound_tightened - new_bound = ub - Mub + 1 - @debug "found global LB tightening $lb -> $new_bound" - if haskey(tree.root.global_tightenings.lower_bounds, j) - if tree.root.global_tightenings.lower_bounds[j] != new_bound - num_tightenings += 1 - end - else - num_tightenings += 1 - end - tree.root.global_tightenings.lower_bounds[j] = new_bound - end - end - node.global_tightenings = num_tightenings - end -end - -""" -Tighten the lower bound using strong convexity of the objective. -""" -function tightening_strong_convexity(tree, x, lower_bound) - μ = tree.root.options[:strong_convexity] - if μ > 0 - @debug "Using strong convexity $μ" - strong_convexity_bound = lower_bound - num_fractional = 0 - for j in tree.root.problem.integer_variables - if x[j] > floor(x[j]) + 1e-6 && x[j] < ceil(x[j]) - 1e-6 - num_fractional += 1 - new_left_increment = μ / 2 * (x[j] - floor(x[j]))^2 - new_right_increment = μ / 2 * (ceil(x[j]) - x[j])^2 - new_increment = min(new_left_increment, new_right_increment) - strong_convexity_bound += new_increment - end - end - @debug "Strong convexity: $lower_bound -> $strong_convexity_bound" - @assert num_fractional == 0 || strong_convexity_bound > lower_bound - lower_bound = strong_convexity_bound - end - return lower_bound -end - - """ Returns the solution vector of the relaxed problem at the node """ diff --git a/src/tightenings.jl b/src/tightenings.jl new file mode 100644 index 000000000..a06777d0e --- /dev/null +++ b/src/tightenings.jl @@ -0,0 +1,247 @@ +""" +Tightening of the bounds at node level. Children node inherit the updated bounds. +""" +function dual_tightening(tree, node, x, dual_gap) + if tree.root.options[:dual_tightening] && isfinite(tree.incumbent) + grad = similar(x) + tree.root.problem.g(grad, x) + num_tightenings = 0 + num_potential_tightenings = 0 + μ = tree.root.options[:strong_convexity] + for j in tree.root.problem.integer_variables + lb_global = get(tree.root.problem.integer_variable_bounds, (j, :greaterthan), -Inf) + ub_global = get(tree.root.problem.integer_variable_bounds, (j, :lessthan), Inf) + lb = get(node.local_bounds.lower_bounds, j, lb_global) + ub = get(node.local_bounds.upper_bounds, j, ub_global) + @assert lb >= lb_global + @assert ub <= ub_global + if lb ≈ ub + # variable already fixed + continue + end + gj = grad[j] + safety_tolerance = 2.0 + rhs = + tree.incumbent - tree.root.problem.f(x) + + safety_tolerance * dual_gap + + sqrt(eps(tree.incumbent)) + if ≈(x[j], lb, atol=tree.options.atol, rtol=tree.options.rtol) + if !isapprox(gj, 0, atol=1e-5) + num_potential_tightenings += 1 + end + if gj > 0 + Mlb = 0 + bound_tightened = true + @debug "starting tightening ub $(rhs)" + while 0.99 * (Mlb * gj + μ / 2 * Mlb^2) <= rhs + Mlb += 1 + if lb + Mlb - 1 == ub + bound_tightened = false + break + end + end + if bound_tightened + new_bound = lb + Mlb - 1 + @debug "found UB tightening $ub -> $new_bound" + node.local_bounds[j, :lessthan] = new_bound + num_tightenings += 1 + if haskey(tree.root.problem.integer_variable_bounds, (j, :lessthan)) + @assert node.local_bounds[j, :lessthan] <= + tree.root.problem.integer_variable_bounds[j, :lessthan] + end + end + end + elseif ≈(x[j], ub, atol=tree.options.atol, rtol=tree.options.rtol) + if !isapprox(gj, 0, atol=1e-5) + num_potential_tightenings += 1 + end + if gj < 0 + Mub = 0 + bound_tightened = true + @debug "starting tightening lb $(rhs)" + while -0.99 * (Mub * gj + μ / 2 * Mub^2) <= rhs + Mub += 1 + if ub - Mub + 1 == lb + bound_tightened = false + break + end + end + if bound_tightened + new_bound = ub - Mub + 1 + @debug "found LB tightening $lb -> $new_bound" + node.local_bounds[j, :greaterthan] = new_bound + num_tightenings += 1 + if haskey(tree.root.problem.integer_variable_bounds, (j, :greaterthan)) + @assert node.local_bounds[j, :greaterthan] >= + tree.root.problem.integer_variable_bounds[j, :greaterthan] + end + end + end + end + end + @debug "# tightenings $num_tightenings" + node.local_tightenings = num_tightenings + node.local_potential_tightenings = num_potential_tightenings + end +end + +""" +Save the gradient of the root solution (i.e. the relaxed solution) and the +corresponding lower and upper bounds. +""" +function store_data_global_tightening(tree, node, x, dual_gap) + if tree.root.options[:global_dual_tightening] && node.std.id == 1 + @debug "storing root node info for tightening" + grad = similar(x) + tree.root.problem.g(grad, x) + safety_tolerance = 2.0 + tree.root.global_tightening_rhs[] = -tree.root.problem.f(x) + safety_tolerance * dual_gap + for j in tree.root.problem.integer_variables + if haskey(tree.root.problem.integer_variable_bounds.upper_bounds, j) + ub = tree.root.problem.integer_variable_bounds[j, :lessthan] + if ≈(x[j], ub, atol=tree.options.atol, rtol=tree.options.rtol) && grad[j] < 0 + tree.root.global_tightening_root_info.upper_bounds[j] = (grad[j], ub) + end + end + if haskey(tree.root.problem.integer_variable_bounds.lower_bounds, j) + lb = tree.root.problem.integer_variable_bounds[j, :greaterthan] + if ≈(x[j], lb, atol=tree.options.atol, rtol=tree.options.rtol) && grad[j] > 0 + tree.root.global_tightening_root_info.lower_bounds[j] = (grad[j], lb) + end + end + end + end +end + +""" +Use the gradient of the root node to tighten the global bounds. +""" +function global_tightening(tree, node) + # new incumbent: check global fixings + if tree.root.options[:global_dual_tightening] && tree.root.updated_incumbent[] + num_tightenings = 0 + rhs = tree.incumbent + tree.root.global_tightening_rhs[] + @assert isfinite(rhs) + for (j, (gj, lb)) in tree.root.global_tightening_root_info.lower_bounds + ub = get(tree.root.problem.integer_variable_bounds.upper_bounds, j, Inf) + ub_new = get(tree.root.global_tightening_root_info.upper_bounds, j, Inf) + ub = min(ub, ub_new) + Mlb = 0 + bound_tightened = true + lb = lb + while Mlb * gj <= rhs + Mlb += 1 + if lb + Mlb - 1 == ub + bound_tightened = false + break + end + end + if bound_tightened + new_bound = lb + Mlb - 1 + @debug "found global UB tightening $ub -> $new_bound" + if haskey(tree.root.global_tightenings.upper_bounds, j) + if tree.root.global_tightenings.upper_bounds[j] != new_bound + num_tightenings += 1 + end + else + num_tightenings += 1 + end + tree.root.global_tightenings.upper_bounds[j] = new_bound + end + end + for (j, (gj, ub)) in tree.root.global_tightening_root_info.upper_bounds + lb = get(tree.root.problem.integer_variable_bounds.lower_bounds, j, -Inf) + lb_new = get(tree.root.global_tightening_root_info.lower_bounds, j, -Inf) + lb = max(lb, lb_new) + Mub = 0 + bound_tightened = true + ub = ub + while -Mub * gj <= rhs + Mub += 1 + if ub - Mub + 1 == lb + bound_tightened = false + break + end + end + if bound_tightened + new_bound = ub - Mub + 1 + @debug "found global LB tightening $lb -> $new_bound" + if haskey(tree.root.global_tightenings.lower_bounds, j) + if tree.root.global_tightenings.lower_bounds[j] != new_bound + num_tightenings += 1 + end + else + num_tightenings += 1 + end + tree.root.global_tightenings.lower_bounds[j] = new_bound + end + end + node.global_tightenings = num_tightenings + end +end + +""" +Tighten the lower bound using strong convexity of the objective. +""" +function tightening_strong_convexity(tree, x, lower_bound) + μ = tree.root.options[:strong_convexity] + if μ > 0 + @debug "Using strong convexity $μ" + strong_convexity_bound = lower_bound + num_fractional = 0 + for j in tree.root.problem.integer_variables + if x[j] > floor(x[j]) + 1e-6 && x[j] < ceil(x[j]) - 1e-6 + num_fractional += 1 + new_left_increment = μ / 2 * (x[j] - floor(x[j]))^2 + new_right_increment = μ / 2 * (ceil(x[j]) - x[j])^2 + new_increment = min(new_left_increment, new_right_increment) + strong_convexity_bound += new_increment + end + end + @debug "Strong convexity: $lower_bound -> $strong_convexity_bound" + @assert num_fractional == 0 || strong_convexity_bound > lower_bound + lower_bound = strong_convexity_bound + end + return lower_bound +end + +""" +Use strong convexity to potentially remove one of the children nodes +""" +function prune_children(tree, node, lower_bound_base, x, vidx) + prune_left = false + prune_right = false + + μ = tree.root.options[:strong_convexity] + if μ > 0 + @debug "Using strong convexity $μ" + for j in tree.root.problem.integer_variables + if vidx == j + continue + end + lower_bound_base += μ / 2 * min((x[j] - floor(x[j]))^2, (ceil(x[j]) - x[j])^2) + 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 + @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 + @debug "prune right, from $(node.lb) -> $new_bound_right, ub $(tree.incumbent), lb $(node.lb)" + prune_right = true + end + @assert !( + (new_bound_left > tree.incumbent + tree.root.options[:dual_gap]) && + (new_bound_right > tree.incumbent + tree.root.options[:dual_gap]) + ) "both sides should not be pruned" + end + + # If both nodes are pruned, when one of them has to be equal to the incumbent. + # Thus, we have proof of optimality by strong convexity. + if prune_left && prune_right + tree.lb = min(new_bound_left, new_bound_right) + end + + return prune_left, prune_right +end From 4b55f60d2bc96f9da134a52306a6eaf6fc88e4ec Mon Sep 17 00:00:00 2001 From: Hendrych Date: Fri, 30 Aug 2024 16:21:47 +0200 Subject: [PATCH 02/45] Tighten bounds using sharpness. --- src/node.jl | 2 +- src/tightenings.jl | 154 ++++++++++++++++++++++++++++++++++++--------- 2 files changed, 126 insertions(+), 30 deletions(-) diff --git a/src/node.jl b/src/node.jl index 66fe189e1..ae60e8a93 100644 --- a/src/node.jl +++ b/src/node.jl @@ -267,7 +267,7 @@ function Bonobo.evaluate_node!(tree::Bonobo.BnBTree, node::FrankWolfeNode) lower_bound = primal - dual_gap # improvement of the lower bound using strong convexity - lower_bound = tightening_strong_convexity(tree, x, lower_bound) + lower_bound = tightening_lowerbound(tree, node, x, lower_bound) # Found an upper bound if is_integer_feasible(tree, x) diff --git a/src/tightenings.jl b/src/tightenings.jl index a06777d0e..09d55e6c4 100644 --- a/src/tightenings.jl +++ b/src/tightenings.jl @@ -181,62 +181,158 @@ function global_tightening(tree, node) end """ -Tighten the lower bound using strong convexity of the objective. +Tighten the lower bound using strong convexity and/or sharpness of the objective. """ -function tightening_strong_convexity(tree, x, lower_bound) +function tightening_lowerbound(tree, node, x, lower_bound) μ = tree.root.options[:strong_convexity] - if μ > 0 - @debug "Using strong convexity $μ" - strong_convexity_bound = lower_bound + M = tree.root.options[:sharpness_constant] + θ = tree.root.options[:sharpness_exponent] + + if μ > 0 || (M > 0 && θ != Inf) + @debug "Tightening lower bound using strong convexity $μ and/or sharpness ($θ, $M)" num_fractional = 0 + bound_improvement = 0.0 for j in tree.root.problem.integer_variables if x[j] > floor(x[j]) + 1e-6 && x[j] < ceil(x[j]) - 1e-6 num_fractional += 1 - new_left_increment = μ / 2 * (x[j] - floor(x[j]))^2 - new_right_increment = μ / 2 * (ceil(x[j]) - x[j])^2 + new_left_increment = (x[j] - floor(x[j]))^2 + new_right_increment = (ceil(x[j]) - x[j])^2 new_increment = min(new_left_increment, new_right_increment) - strong_convexity_bound += new_increment + bound_improvement += new_increment end end - @debug "Strong convexity: $lower_bound -> $strong_convexity_bound" - @assert num_fractional == 0 || strong_convexity_bound > lower_bound - lower_bound = strong_convexity_bound + strong_convexity_bound = -Inf + sharpness_bound = -Inf + + # strong convexity + if μ > 0 + @debug "Using strong convexity $μ" + strong_convexity = lower_bound + μ / 2 * bound_improvement + @debug "Strong convexity: $lower_bound -> $strong_convexity_bound" + @assert num_fractional == 0 || strong_convexity_bound > lower_bound + end + + # sharpness + if M > 0 && θ != Inf + + @debug "Using sharpness θ=$θ and M=$M" + fx = tree.root.problem.f(x) + + # build the root LMO + free_model(tree.root.problem.tlmo.blmo) + build_LMO( + tree.root.problem.tlmo, + tree.root.problem.integer_variable_bounds, + IntegerBounds(), + tree.root.problem.integer_variables, + ) + + gradient = zeros(tree.root.problem.nvars) + tree.root.problem.g(gradient, x) + v = compute_extreme_point(tree.root.problem.tlmo, gradient) + root_dual_gap = FrankWolfe.fast_dot(gradient, x) - FrankWolfe.fast_dot(gradient, v) + + free_model(tree.root.problem.tlmo.blmo) + build_LMO( + tree.root.problem.tlmo, + tree.root.problem.integer_variable_bounds, + node.local_bounds, + tree.root.problem.integer_variables, + ) + + sharpness_bound = M^(- 1 / θ) * (sqrt(bound_improvement) - M * root_dual_gap^θ)^(1 / θ) + fx - root_dual_gap + @debug "Sharpness: $lower_bound -> $sharpness_bound" + @assert num_fractional == 0 || sharpness_bound > lower_bound + end + + lower_bound = max(strong_convexity_bound, sharpness_bound) end return lower_bound end """ -Use strong convexity to potentially remove one of the children nodes +Use strong convexity and/or sharpness to potentially remove one of the children nodes """ function prune_children(tree, node, lower_bound_base, x, vidx) prune_left = false prune_right = false μ = tree.root.options[:strong_convexity] - if μ > 0 - @debug "Using strong convexity $μ" + M = tree.root.options[:sharpness_constant] + θ = tree.root.options[:sharpness_exponent] + + strong_convexity_bound + + if μ > 0 || (M > 0 && θ != Inf) + bound_improvement = 0.0 for j in tree.root.problem.integer_variables if vidx == j continue end - lower_bound_base += μ / 2 * min((x[j] - floor(x[j]))^2, (ceil(x[j]) - x[j])^2) - 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 - @debug "prune left, from $(node.lb) -> $new_bound_left, ub $(tree.incumbent), lb $(node.lb)" - prune_left = true + bound_improvement += min((x[j] - floor(x[j]))^2, (ceil(x[j]) - x[j])^2) end - if new_bound_right > tree.incumbent - @debug "prune right, from $(node.lb) -> $new_bound_right, ub $(tree.incumbent), lb $(node.lb)" - prune_right = true + new_bound_left = bound_improvement + (x[vidx] - floor(x[vidx]))^2 + new_bound_right = bound_improvement + (ceil(x[vidx]) - x[vidx])^2 + + # strong convexity + if μ > 0 + new_bound_left = lower_bound_base + μ / 2 * new_bound_left + new_bound_right = lower_bound_base + μ / 2 * new_bound_right + + 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 + @debug "prune right, from $(node.lb) -> $new_bound_right, ub $(tree.incumbent), lb $(node.lb)" + prune_right = true + end + @assert !( + (new_bound_left > tree.incumbent + tree.root.options[:dual_gap]) && + (new_bound_right > tree.incumbent + tree.root.options[:dual_gap]) + ) + # sharpness + elseif M > 0 && θ != Inf + x = tree.root.problem.f(x) + # build the root LMO + free_model(tree.root.problem.tlmo.blmo) + build_LMO( + tree.root.problem.tlmo, + tree.root.problem.integer_variable_bounds, + IntegerBounds(), + tree.root.problem.integer_variables, + ) + gradient = zeros(tree.root.problem.nvars) + tree.root.problem.g(gradient, x) + v = compute_extreme_point(tree.root.problem.tlmo, gradient) + root_dual_gap = FrankWolfe.fast_dot(gradient, x) - FrankWolfe.fast_dot(gradient, v) + + free_model(tree.root.problem.tlmo.blmo) + build_LMO( + tree.root.problem.tlmo, + tree.root.problem.integer_variable_bounds, + node.local_bounds, + tree.root.problem.integer_variables, + ) + + new_bound_left = M^(- 1 / θ) * (sqrt(new_bound_left) - M * root_dual_gap^θ)^(1 / θ) + fx - root_dual_gap + new_bound_right = M^(- 1 / θ) * (sqrt(new_bound_right) - M * root_dual_gap^θ)^(1 / θ) + fx - root_dual_gap + + 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 + @debug "prune right, from $(node.lb) -> $new_bound_right, ub $(tree.incumbent), lb $(node.lb)" + prune_right = true + end + @assert !( + (new_bound_left > tree.incumbent + tree.root.options[:dual_gap]) && + (new_bound_right > tree.incumbent + tree.root.options[:dual_gap]) + ) end - @assert !( - (new_bound_left > tree.incumbent + tree.root.options[:dual_gap]) && - (new_bound_right > tree.incumbent + tree.root.options[:dual_gap]) - ) "both sides should not be pruned" end - + # If both nodes are pruned, when one of them has to be equal to the incumbent. # Thus, we have proof of optimality by strong convexity. if prune_left && prune_right From 7f9b4c73a5b20fe7806e7ac946f3e28c7cef6976 Mon Sep 17 00:00:00 2001 From: Hendrych Date: Mon, 2 Sep 2024 12:25:30 +0200 Subject: [PATCH 03/45] Add sharpness constant and exponent to the interface. --- src/MOI_bounded_oracle.jl | 4 ++++ src/interface.jl | 6 ++++++ src/managed_blmo.jl | 4 ++++ 3 files changed, 14 insertions(+) diff --git a/src/MOI_bounded_oracle.jl b/src/MOI_bounded_oracle.jl index f3e7f74fc..d72c5786b 100644 --- a/src/MOI_bounded_oracle.jl +++ b/src/MOI_bounded_oracle.jl @@ -629,6 +629,8 @@ function solve( global_dual_tightening=true, bnb_callback=nothing, strong_convexity=0.0, + sharpness_constant = 0.0, + sharpness_exponent = Inf, domain_oracle=x -> true, start_solution=nothing, fw_verbose=false, @@ -664,6 +666,8 @@ function solve( global_dual_tightening=global_dual_tightening, bnb_callback=bnb_callback, strong_convexity=strong_convexity, + sharpness_constant=sharpness_constant, + sharpness_exponent=sharpness_exponent, domain_oracle=domain_oracle, start_solution=start_solution, fw_verbose=fw_verbose, diff --git a/src/interface.jl b/src/interface.jl index 65fe63ad9..e3a29ff90 100644 --- a/src/interface.jl +++ b/src/interface.jl @@ -43,6 +43,8 @@ dual_tightening - whether to use dual tightening techniques (make sure yo 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 +sharpness_constant - the constant M > 0 for (θ, M)-sharpness. +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 on the Line Search method, you might have to provide the domain oracle to it, too. @@ -83,6 +85,8 @@ function solve( global_dual_tightening=true, bnb_callback=nothing, strong_convexity=0.0, + sharpness_constant = 0.0, + sharpness_exponent = Inf, domain_oracle=x -> true, start_solution=nothing, fw_verbose=false, @@ -196,6 +200,8 @@ function solve( :max_fw_iter => max_fw_iter, :print_iter => print_iter, :strong_convexity => strong_convexity, + :sharpness_constant => sharpness_constant, + :sharpness_exponent => sharpness_exponent, :time_limit => time_limit, :usePostsolve => use_postsolve, :variant => variant, diff --git a/src/managed_blmo.jl b/src/managed_blmo.jl index e2d5e7ebe..21d6ae7ac 100644 --- a/src/managed_blmo.jl +++ b/src/managed_blmo.jl @@ -269,6 +269,8 @@ function solve( global_dual_tightening=true, bnb_callback=nothing, strong_convexity=0.0, + sharpness_constant = 0.0, + sharpness_exponent = Inf, domain_oracle=x -> true, start_solution=nothing, fw_verbose=false, @@ -306,6 +308,8 @@ function solve( global_dual_tightening=global_dual_tightening, bnb_callback=bnb_callback, strong_convexity=strong_convexity, + sharpness_constant=sharpness_constant, + sharpness_exponent=sharpness_exponent, domain_oracle=domain_oracle, start_solution=start_solution, fw_verbose=fw_verbose, From 5819c747910eb7fdc64d2f54d3ce87f3934fb22f Mon Sep 17 00:00:00 2001 From: Hendrych Date: Mon, 2 Sep 2024 13:52:46 +0200 Subject: [PATCH 04/45] Minor fixes. --- src/tightenings.jl | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/src/tightenings.jl b/src/tightenings.jl index 09d55e6c4..089c76401 100644 --- a/src/tightenings.jl +++ b/src/tightenings.jl @@ -201,15 +201,15 @@ function tightening_lowerbound(tree, node, x, lower_bound) bound_improvement += new_increment end end - strong_convexity_bound = -Inf + strong_convexity_bound = lower_bound sharpness_bound = -Inf # strong convexity if μ > 0 @debug "Using strong convexity $μ" - strong_convexity = lower_bound + μ / 2 * bound_improvement + strong_convexity_bound += μ / 2 * bound_improvement @debug "Strong convexity: $lower_bound -> $strong_convexity_bound" - @assert num_fractional == 0 || strong_convexity_bound > lower_bound + @assert num_fractional == 0 || strong_convexity_bound > lower_bound end # sharpness @@ -261,8 +261,6 @@ function prune_children(tree, node, lower_bound_base, x, vidx) M = tree.root.options[:sharpness_constant] θ = tree.root.options[:sharpness_exponent] - strong_convexity_bound - if μ > 0 || (M > 0 && θ != Inf) bound_improvement = 0.0 for j in tree.root.problem.integer_variables From 81000b6c4db4d4149b973ff15ce0322d1865ce55 Mon Sep 17 00:00:00 2001 From: Hendrych Date: Mon, 2 Sep 2024 19:00:51 +0200 Subject: [PATCH 05/45] Tests for strong convexity and sharpness on the Normbox example. --- test/interface_test.jl | 57 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 57 insertions(+) diff --git a/test/interface_test.jl b/test/interface_test.jl index 572f0b557..cbd812c13 100644 --- a/test/interface_test.jl +++ b/test/interface_test.jl @@ -117,6 +117,63 @@ end end end +@testset "Normbox - strong convexity and sharpness" begin + 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 + + @testset "Strong convexity" begin + int_vars = collect(1:n) + lbs = zeros(n) + ubs = ones(n) + + sblmo = Boscia.CubeSimpleBLMO(lbs, ubs, int_vars) + μ = 1.0 + + x, _, result = + Boscia.solve(f, grad!, sblmo, lbs[int_vars], ubs[int_vars], int_vars, n, strong_convexity=μ) + + @test x == round.(diffi) + @test isapprox(f(x), f(result[:raw_solution]), atol=1e-6, rtol=1e-3) + end + + @testset "Sharpness" begin + int_vars = collect(1:n) + lbs = zeros(n) + ubs = ones(n) + + sblmo = Boscia.CubeSimpleBLMO(lbs, ubs, int_vars) + θ = 1/2 + M = 2.0 + + x, _, result = + Boscia.solve(f, grad!, sblmo, lbs[int_vars], ubs[int_vars], int_vars, n, sharpness_constant=M, sharpness_exponent=θ) + + @test x == round.(diffi) + @test isapprox(f(x), f(result[:raw_solution]), atol=1e-6, rtol=1e-3) + end + + @testset "Strong convexity and sharpness" begin + int_vars = collect(1:n) + lbs = zeros(n) + ubs = ones(n) + + sblmo = Boscia.CubeSimpleBLMO(lbs, ubs, int_vars) + μ = 1.0 + θ = 1/2 + M = 2.0 + + x, _, result = + Boscia.solve(f, grad!, sblmo, lbs[int_vars], ubs[int_vars], int_vars, n, strong_convexity=μ, sharpness_constant=M, sharpness_exponent=θ) + + @test x == round.(diffi) + @test isapprox(f(x), f(result[:raw_solution]), atol=1e-6, rtol=1e-3) + end +end + @testset "Start with Active Set" begin function f(x) From df136693ea8c02a8b2e3870bfaee74820928cc4e Mon Sep 17 00:00:00 2001 From: Hendrych Date: Mon, 2 Sep 2024 19:09:51 +0200 Subject: [PATCH 06/45] Use the node dual gap instead of root dual gap. --- src/tightenings.jl | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/src/tightenings.jl b/src/tightenings.jl index 089c76401..ddb01c8ce 100644 --- a/src/tightenings.jl +++ b/src/tightenings.jl @@ -240,7 +240,8 @@ function tightening_lowerbound(tree, node, x, lower_bound) tree.root.problem.integer_variables, ) - sharpness_bound = M^(- 1 / θ) * (sqrt(bound_improvement) - M * root_dual_gap^θ)^(1 / θ) + fx - root_dual_gap + #sharpness_bound = M^(- 1 / θ) * (sqrt(bound_improvement) - M * root_dual_gap^θ)^(1 / θ) + fx - root_dual_gap + sharpness_bound = M^(- 1 / θ) * (sqrt(bound_improvement) - M * node.dual_gap^θ)^(1 / θ) + fx - node.dual_gap @debug "Sharpness: $lower_bound -> $sharpness_bound" @assert num_fractional == 0 || sharpness_bound > lower_bound end @@ -291,7 +292,7 @@ function prune_children(tree, node, lower_bound_base, x, vidx) ) # sharpness elseif M > 0 && θ != Inf - x = tree.root.problem.f(x) + fx = tree.root.problem.f(x) # build the root LMO free_model(tree.root.problem.tlmo.blmo) build_LMO( @@ -313,8 +314,11 @@ function prune_children(tree, node, lower_bound_base, x, vidx) tree.root.problem.integer_variables, ) - new_bound_left = M^(- 1 / θ) * (sqrt(new_bound_left) - M * root_dual_gap^θ)^(1 / θ) + fx - root_dual_gap - new_bound_right = M^(- 1 / θ) * (sqrt(new_bound_right) - M * root_dual_gap^θ)^(1 / θ) + fx - root_dual_gap + #new_bound_left = M^(- 1 / θ) * (sqrt(new_bound_left) - M * root_dual_gap^θ)^(1 / θ) + fx - root_dual_gap + #new_bound_right = M^(- 1 / θ) * (sqrt(new_bound_right) - M * root_dual_gap^θ)^(1 / θ) + fx - root_dual_gap + + new_bound_left = M^(- 1 / θ) * (sqrt(new_bound_left) - M * node.dual_gap^θ)^(1 / θ) + fx - node.dual_gap + new_bound_right = M^(- 1 / θ) * (sqrt(new_bound_right) - M * node.dual_gap^θ)^(1 / θ) + fx - node.dual_gap if new_bound_left > tree.incumbent @debug "prune left, from $(node.lb) -> $new_bound_left, ub $(tree.incumbent), lb $(node.lb)" From a17f99362c32b0c943c77d7426c527fe5a367b77 Mon Sep 17 00:00:00 2001 From: Hendrych Date: Tue, 3 Sep 2024 12:46:34 +0200 Subject: [PATCH 07/45] Make sure that we collect the statistics, even if in case of timeout. --- src/custom_bonobo.jl | 19 +++++++++++++++++++ src/interface.jl | 8 ++++---- 2 files changed, 23 insertions(+), 4 deletions(-) diff --git a/src/custom_bonobo.jl b/src/custom_bonobo.jl index 0fcf91775..46780811e 100644 --- a/src/custom_bonobo.jl +++ b/src/custom_bonobo.jl @@ -75,6 +75,25 @@ function Bonobo.optimize!( Bonobo.branch!(tree, node) callback(tree, node) end + # To make sure that we collect the statistics in case the time limit is crashed. + if haskey(tree.root.result, :global_tightenings) + y = Bonobo.get_solution(tree) + vertex_storage = FrankWolfe.DeletedVertexStorage(typeof(y)[], 1) + dummy_node = FrankWolfeNode( + NodeInfo(-1, Inf, Inf), + FrankWolfe.ActiveSet([(1.0, y)]), + vertex_storage, + IntegerBounds(), + 1, + 1e-3, + Millisecond(0), + 0, + 0, + 0, + 0.0, + ) + callback(tree, dummy_node, node_infeasible=true) + end return Bonobo.sort_solutions!(tree.solutions, tree.sense) end diff --git a/src/interface.jl b/src/interface.jl index e3a29ff90..24f6be846 100644 --- a/src/interface.jl +++ b/src/interface.jl @@ -209,6 +209,7 @@ function solve( :heuristics => heuristics, :heu_ncalls => 0, ), + result=Dict{Symbol,Any}(), ), branch_strategy=branching_strategy, dual_gap_limit=rel_dual_gap, @@ -258,7 +259,6 @@ function solve( list_discarded_set_size_cb = Int[] fw_iterations = Int[] node_level = Int[] - result = Dict{Symbol,Any}() lmo_calls_per_layer = Vector{Vector{Int}}() active_set_size_per_layer = Vector{Vector{Int}}() discarded_set_size_per_layer = Vector{Vector{Int}}() @@ -279,7 +279,7 @@ function solve( fw_iterations, list_active_set_size_cb, list_discarded_set_size_cb, - result, + tree.root.result, lmo_calls_per_layer, active_set_size_per_layer, discarded_set_size_per_layer, @@ -299,7 +299,7 @@ function solve( Bonobo.optimize!(tree; callback=bnb_callback) - x = postsolve(tree, result, time_ref, verbose, max_iteration_post) + x = postsolve(tree, tree.root.result, time_ref, verbose, max_iteration_post) # Check solution and polish x_polished = x @@ -322,7 +322,7 @@ function solve( end println() # cleaner output - return x, tree.root.problem.tlmo, result + return x, tree.root.problem.tlmo, tree.root.result end """ From 06a0e25eaba8a718928078c8083e4f746bd8c1f5 Mon Sep 17 00:00:00 2001 From: Hendrych Date: Tue, 3 Sep 2024 12:51:36 +0200 Subject: [PATCH 08/45] We only need this if the key is not in the result dictionary. --- src/custom_bonobo.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/custom_bonobo.jl b/src/custom_bonobo.jl index 46780811e..d86094e49 100644 --- a/src/custom_bonobo.jl +++ b/src/custom_bonobo.jl @@ -76,7 +76,7 @@ function Bonobo.optimize!( callback(tree, node) end # To make sure that we collect the statistics in case the time limit is crashed. - if haskey(tree.root.result, :global_tightenings) + if !haskey(tree.root.result, :global_tightenings) y = Bonobo.get_solution(tree) vertex_storage = FrankWolfe.DeletedVertexStorage(typeof(y)[], 1) dummy_node = FrankWolfeNode( From c8add285fe5bbc9c6d4c58bc299110228569d71d Mon Sep 17 00:00:00 2001 From: Hendrych Date: Tue, 3 Sep 2024 12:56:08 +0200 Subject: [PATCH 09/45] Move the BnB Callback into the appriopiate file. --- src/callbacks.jl | 260 ++++++++++++++++++++++++++++++++++++++++++++++- src/interface.jl | 248 -------------------------------------------- 2 files changed, 259 insertions(+), 249 deletions(-) diff --git a/src/callbacks.jl b/src/callbacks.jl index dcb14f3ee..8da2eff98 100644 --- a/src/callbacks.jl +++ b/src/callbacks.jl @@ -1,4 +1,11 @@ -# FW callback +""" +Frank-Wolfe Callback. + +Is called in every Frank-Wolfe iteration. +Node evaluation can be dymamically stopped here. +Time limit is checked. +If the vertex is providing a better incumbent, it is added as solution. +""" function build_FW_callback( tree, min_number_lower, @@ -94,3 +101,254 @@ function build_FW_callback( return true end end + +""" +Branch-and-Bound Callback. +Collects statistics and prints logs if verbose is turned on. + +Output of Boscia + + iter : current iteration of Boscia + node id : current node id + lower bound : tree_lb(tree) + incumbent : tree.incumbent + gap : tree.incumbent-tree_lb(tree) + rel. gap : dual_gap/tree.incumbent + time : total time of Boscia + time/nodes : average time per node + FW time : time spent in FW + LMO time : time used by LMO + LMO calls : number of compute_extreme_point calls in FW + FW iterations : number of iterations in FW +""" +function build_bnb_callback( + tree, + time_ref, + list_lb_cb, + list_ub_cb, + list_time_cb, + list_num_nodes_cb, + list_lmo_calls_cb, + verbose, + fw_iterations, + list_active_set_size_cb, + list_discarded_set_size_cb, + result, + lmo_calls_per_layer, + active_set_size_per_layer, + discarded_set_size_per_layer, + node_level, + baseline_callback, + local_tightenings, + global_tightenings, + local_potential_tightenings, + num_bin, + num_int, +) + iteration = 0 + + headers = [ + " ", + "Iter", + "Open", + "Bound", + "Incumbent", + "Gap (abs)", + "Gap (rel)", + "Time (s)", + "Nodes/sec", + "FW (ms)", + "LMO (ms)", + "LMO (calls c)", + "FW (its)", + "#activeset", + "#shadow", + ] + format_string = "%1s %5i %5i %14e %14e %14e %14e %14e %14e %12i %10i %14i %10i %8i %8i\n" + print_iter = get(tree.root.options, :print_iter, 100) + + if verbose + FrankWolfe.print_callback(headers, format_string, print_header=true) + end + return function callback( + tree, + node; + worse_than_incumbent=false, + node_infeasible=false, + lb_update=false, + ) + if baseline_callback !== nothing + baseline_callback( + tree, + node, + worse_than_incumbent=worse_than_incumbent, + node_infeasible=node_infeasible, + lb_update=lb_update, + ) + end + if !node_infeasible + #update lower bound + if lb_update == true + tree.node_queue[node.id] = (node.lb, node.id) + _, prio = peek(tree.node_queue) + @assert tree.lb <= prio[1] + tree.lb = min(minimum([prio[2][1] for prio in tree.node_queue]), tree.incumbent) + end + push!(list_ub_cb, tree.incumbent) + push!(list_num_nodes_cb, tree.num_nodes) + push!(node_level, node.level) + iteration += 1 + if tree.lb == -Inf && isempty(tree.nodes) + tree.lb = node.lb + end + + time = float(Dates.value(Dates.now() - time_ref)) + push!(list_time_cb, time) + + if tree.root.options[:time_limit] < Inf + if time / 1000.0 ≥ tree.root.options[:time_limit] + @assert tree.root.problem.solving_stage == SOLVING + tree.root.problem.solving_stage = TIME_LIMIT_REACHED + end + end + + fw_time = Dates.value(node.fw_time) + fw_iter = if !isempty(fw_iterations) + fw_iterations[end] + else + 0 + end + if !isempty(tree.root.problem.tlmo.optimizing_times) + LMO_time = sum(1000 * tree.root.problem.tlmo.optimizing_times) + empty!(tree.root.problem.tlmo.optimizing_times) + else + LMO_time = 0 + end + LMO_calls_c = tree.root.problem.tlmo.ncalls + push!(list_lmo_calls_cb, copy(LMO_calls_c)) + + if !isempty(tree.node_queue) + p_lb = tree.lb + tree.lb = min(minimum([prio[2][1] for prio in tree.node_queue]), tree.incumbent) + @assert p_lb <= tree.lb + tree.root.options[:dual_gap] "p_lb <= tree.lb + tree.root.options[:dual_gap] $(p_lb) <= $(tree.lb + tree.root.options[:dual_gap])" + end + # correct lower bound if necessary + tree.lb = tree_lb(tree) + dual_gap = tree.incumbent - tree_lb(tree) + push!(list_lb_cb, tree_lb(tree)) + active_set_size = length(node.active_set) + discarded_set_size = length(node.discarded_vertices.storage) + push!(list_active_set_size_cb, active_set_size) + push!(list_discarded_set_size_cb, discarded_set_size) + nodes_left = length(tree.nodes) + if tree.root.updated_incumbent[] + incumbent_updated = "*" + else + incumbent_updated = " " + end + if verbose && ( + mod(iteration, print_iter) == 0 || + iteration == 1 || + Bonobo.terminated(tree) || + tree.root.updated_incumbent[] + ) + if (mod(iteration, print_iter * 40) == 0) + FrankWolfe.print_callback(headers, format_string, print_header=true) + end + FrankWolfe.print_callback( + ( + incumbent_updated, + iteration, + nodes_left, + tree_lb(tree), + tree.incumbent, + dual_gap, + relative_gap(tree.incumbent, tree_lb(tree)), + time / 1000.0, + tree.num_nodes / time * 1000.0, + fw_time, + LMO_time, + tree.root.problem.tlmo.ncalls, + fw_iter, + active_set_size, + discarded_set_size, + ), + format_string, + print_header=false, + ) + tree.root.updated_incumbent[] = false + end + # lmo calls per layer + if length(list_lmo_calls_cb) > 1 + LMO_calls = list_lmo_calls_cb[end] - list_lmo_calls_cb[end-1] + else + LMO_calls = list_lmo_calls_cb[end] + end + if length(lmo_calls_per_layer) < node.level + push!(lmo_calls_per_layer, [LMO_calls]) + push!(active_set_size_per_layer, [active_set_size]) + push!(discarded_set_size_per_layer, [discarded_set_size]) + else + push!(lmo_calls_per_layer[node.level], LMO_calls) + push!(active_set_size_per_layer[node.level], active_set_size) + push!(discarded_set_size_per_layer[node.level], discarded_set_size) + end + + # add tightenings + push!(global_tightenings, node.global_tightenings) + push!(local_tightenings, node.local_tightenings) + push!(local_potential_tightenings, node.local_potential_tightenings) + @assert node.local_potential_tightenings <= num_bin + num_int + @assert node.local_tightenings <= num_bin + num_int + @assert node.global_tightenings <= num_bin + num_int + end + # update current_node_id + if !Bonobo.terminated(tree) + tree.root.current_node_id[] = + Bonobo.get_next_node(tree, tree.options.traverse_strategy).id + end + + if Bonobo.terminated(tree) + Bonobo.sort_solutions!(tree.solutions, tree.sense) + x = Bonobo.get_solution(tree) + # 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 + + # If the tree is empty, incumbent and solution should be the same! + if isempty(tree.nodes) + @assert isapprox(tree.incumbent, primal_value) + end + + result[:number_nodes] = tree.num_nodes + result[:lmo_calls] = tree.root.problem.tlmo.ncalls + result[:heu_lmo_calls] = tree.root.options[:heu_ncalls] + result[:list_num_nodes] = list_num_nodes_cb + result[:list_lmo_calls_acc] = list_lmo_calls_cb + result[:list_active_set_size] = list_active_set_size_cb + result[:list_discarded_set_size] = list_discarded_set_size_cb + result[:list_lb] = list_lb_cb + result[:list_ub] = list_ub_cb + result[:list_time] = list_time_cb + result[:lmo_calls_per_layer] = lmo_calls_per_layer + result[:active_set_size_per_layer] = active_set_size_per_layer + result[:discarded_set_size_per_layer] = discarded_set_size_per_layer + result[:node_level] = node_level + result[:global_tightenings] = global_tightenings + result[:local_tightenings] = local_tightenings + result[:local_potential_tightenings] = local_potential_tightenings + + if verbose + FrankWolfe.print_callback(headers, format_string, print_footer=true) + println() + end + end + end +end diff --git a/src/interface.jl b/src/interface.jl index 24f6be846..920ad881f 100644 --- a/src/interface.jl +++ b/src/interface.jl @@ -325,254 +325,6 @@ function solve( return x, tree.root.problem.tlmo, tree.root.result end -""" -Output of Boscia - - iter : current iteration of Boscia - node id : current node id - lower bound : tree_lb(tree) - incumbent : tree.incumbent - gap : tree.incumbent-tree_lb(tree) - rel. gap : dual_gap/tree.incumbent - time : total time of Boscia - time/nodes : average time per node - FW time : time spent in FW - LMO time : time used by LMO - LMO calls : number of compute_extreme_point calls in FW - FW iterations : number of iterations in FW -""" -function build_bnb_callback( - tree, - time_ref, - list_lb_cb, - list_ub_cb, - list_time_cb, - list_num_nodes_cb, - list_lmo_calls_cb, - verbose, - fw_iterations, - list_active_set_size_cb, - list_discarded_set_size_cb, - result, - lmo_calls_per_layer, - active_set_size_per_layer, - discarded_set_size_per_layer, - node_level, - baseline_callback, - local_tightenings, - global_tightenings, - local_potential_tightenings, - num_bin, - num_int, -) - iteration = 0 - - headers = [ - " ", - "Iter", - "Open", - "Bound", - "Incumbent", - "Gap (abs)", - "Gap (rel)", - "Time (s)", - "Nodes/sec", - "FW (ms)", - "LMO (ms)", - "LMO (calls c)", - "FW (its)", - "#activeset", - "#shadow", - ] - format_string = "%1s %5i %5i %14e %14e %14e %14e %14e %14e %12i %10i %14i %10i %8i %8i\n" - print_iter = get(tree.root.options, :print_iter, 100) - - if verbose - FrankWolfe.print_callback(headers, format_string, print_header=true) - end - return function callback( - tree, - node; - worse_than_incumbent=false, - node_infeasible=false, - lb_update=false, - ) - if baseline_callback !== nothing - baseline_callback( - tree, - node, - worse_than_incumbent=worse_than_incumbent, - node_infeasible=node_infeasible, - lb_update=lb_update, - ) - end - if !node_infeasible - #update lower bound - if lb_update == true - tree.node_queue[node.id] = (node.lb, node.id) - _, prio = peek(tree.node_queue) - @assert tree.lb <= prio[1] - tree.lb = min(minimum([prio[2][1] for prio in tree.node_queue]), tree.incumbent) - end - push!(list_ub_cb, tree.incumbent) - push!(list_num_nodes_cb, tree.num_nodes) - push!(node_level, node.level) - iteration += 1 - if tree.lb == -Inf && isempty(tree.nodes) - tree.lb = node.lb - end - - time = float(Dates.value(Dates.now() - time_ref)) - push!(list_time_cb, time) - - if tree.root.options[:time_limit] < Inf - if time / 1000.0 ≥ tree.root.options[:time_limit] - @assert tree.root.problem.solving_stage == SOLVING - tree.root.problem.solving_stage = TIME_LIMIT_REACHED - end - end - - fw_time = Dates.value(node.fw_time) - fw_iter = if !isempty(fw_iterations) - fw_iterations[end] - else - 0 - end - if !isempty(tree.root.problem.tlmo.optimizing_times) - LMO_time = sum(1000 * tree.root.problem.tlmo.optimizing_times) - empty!(tree.root.problem.tlmo.optimizing_times) - else - LMO_time = 0 - end - LMO_calls_c = tree.root.problem.tlmo.ncalls - push!(list_lmo_calls_cb, copy(LMO_calls_c)) - - if !isempty(tree.node_queue) - p_lb = tree.lb - tree.lb = min(minimum([prio[2][1] for prio in tree.node_queue]), tree.incumbent) - @assert p_lb <= tree.lb + tree.root.options[:dual_gap] "p_lb <= tree.lb + tree.root.options[:dual_gap] $(p_lb) <= $(tree.lb + tree.root.options[:dual_gap])" - end - # correct lower bound if necessary - tree.lb = tree_lb(tree) - dual_gap = tree.incumbent - tree_lb(tree) - push!(list_lb_cb, tree_lb(tree)) - active_set_size = length(node.active_set) - discarded_set_size = length(node.discarded_vertices.storage) - push!(list_active_set_size_cb, active_set_size) - push!(list_discarded_set_size_cb, discarded_set_size) - nodes_left = length(tree.nodes) - if tree.root.updated_incumbent[] - incumbent_updated = "*" - else - incumbent_updated = " " - end - if verbose && ( - mod(iteration, print_iter) == 0 || - iteration == 1 || - Bonobo.terminated(tree) || - tree.root.updated_incumbent[] - ) - if (mod(iteration, print_iter * 40) == 0) - FrankWolfe.print_callback(headers, format_string, print_header=true) - end - FrankWolfe.print_callback( - ( - incumbent_updated, - iteration, - nodes_left, - tree_lb(tree), - tree.incumbent, - dual_gap, - relative_gap(tree.incumbent, tree_lb(tree)), - time / 1000.0, - tree.num_nodes / time * 1000.0, - fw_time, - LMO_time, - tree.root.problem.tlmo.ncalls, - fw_iter, - active_set_size, - discarded_set_size, - ), - format_string, - print_header=false, - ) - tree.root.updated_incumbent[] = false - end - # lmo calls per layer - if length(list_lmo_calls_cb) > 1 - LMO_calls = list_lmo_calls_cb[end] - list_lmo_calls_cb[end-1] - else - LMO_calls = list_lmo_calls_cb[end] - end - if length(lmo_calls_per_layer) < node.level - push!(lmo_calls_per_layer, [LMO_calls]) - push!(active_set_size_per_layer, [active_set_size]) - push!(discarded_set_size_per_layer, [discarded_set_size]) - else - push!(lmo_calls_per_layer[node.level], LMO_calls) - push!(active_set_size_per_layer[node.level], active_set_size) - push!(discarded_set_size_per_layer[node.level], discarded_set_size) - end - - # add tightenings - push!(global_tightenings, node.global_tightenings) - push!(local_tightenings, node.local_tightenings) - push!(local_potential_tightenings, node.local_potential_tightenings) - @assert node.local_potential_tightenings <= num_bin + num_int - @assert node.local_tightenings <= num_bin + num_int - @assert node.global_tightenings <= num_bin + num_int - end - # update current_node_id - if !Bonobo.terminated(tree) - tree.root.current_node_id[] = - Bonobo.get_next_node(tree, tree.options.traverse_strategy).id - end - - if Bonobo.terminated(tree) - Bonobo.sort_solutions!(tree.solutions, tree.sense) - x = Bonobo.get_solution(tree) - # 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 - - # If the tree is empty, incumbent and solution should be the same! - if isempty(tree.nodes) - @assert isapprox(tree.incumbent, primal_value) - end - - result[:number_nodes] = tree.num_nodes - result[:lmo_calls] = tree.root.problem.tlmo.ncalls - result[:heu_lmo_calls] = tree.root.options[:heu_ncalls] - result[:list_num_nodes] = list_num_nodes_cb - result[:list_lmo_calls_acc] = list_lmo_calls_cb - result[:list_active_set_size] = list_active_set_size_cb - result[:list_discarded_set_size] = list_discarded_set_size_cb - result[:list_lb] = list_lb_cb - result[:list_ub] = list_ub_cb - result[:list_time] = list_time_cb - result[:lmo_calls_per_layer] = lmo_calls_per_layer - result[:active_set_size_per_layer] = active_set_size_per_layer - result[:discarded_set_size_per_layer] = discarded_set_size_per_layer - result[:node_level] = node_level - result[:global_tightenings] = global_tightenings - result[:local_tightenings] = local_tightenings - result[:local_potential_tightenings] = local_potential_tightenings - - if verbose - FrankWolfe.print_callback(headers, format_string, print_footer=true) - println() - end - end - end -end - """ postsolve(tree, result, time_ref, verbose) From 0fad6d282faa90854d7de63c9670471d9ff18261 Mon Sep 17 00:00:00 2001 From: Hendrych Date: Tue, 3 Sep 2024 13:04:58 +0200 Subject: [PATCH 10/45] We don't need the root dual gap. --- src/tightenings.jl | 46 ---------------------------------------------- 1 file changed, 46 deletions(-) diff --git a/src/tightenings.jl b/src/tightenings.jl index ddb01c8ce..2d9ead276 100644 --- a/src/tightenings.jl +++ b/src/tightenings.jl @@ -218,29 +218,6 @@ function tightening_lowerbound(tree, node, x, lower_bound) @debug "Using sharpness θ=$θ and M=$M" fx = tree.root.problem.f(x) - # build the root LMO - free_model(tree.root.problem.tlmo.blmo) - build_LMO( - tree.root.problem.tlmo, - tree.root.problem.integer_variable_bounds, - IntegerBounds(), - tree.root.problem.integer_variables, - ) - - gradient = zeros(tree.root.problem.nvars) - tree.root.problem.g(gradient, x) - v = compute_extreme_point(tree.root.problem.tlmo, gradient) - root_dual_gap = FrankWolfe.fast_dot(gradient, x) - FrankWolfe.fast_dot(gradient, v) - - free_model(tree.root.problem.tlmo.blmo) - build_LMO( - tree.root.problem.tlmo, - tree.root.problem.integer_variable_bounds, - node.local_bounds, - tree.root.problem.integer_variables, - ) - - #sharpness_bound = M^(- 1 / θ) * (sqrt(bound_improvement) - M * root_dual_gap^θ)^(1 / θ) + fx - root_dual_gap sharpness_bound = M^(- 1 / θ) * (sqrt(bound_improvement) - M * node.dual_gap^θ)^(1 / θ) + fx - node.dual_gap @debug "Sharpness: $lower_bound -> $sharpness_bound" @assert num_fractional == 0 || sharpness_bound > lower_bound @@ -293,29 +270,6 @@ function prune_children(tree, node, lower_bound_base, x, vidx) # sharpness elseif M > 0 && θ != Inf fx = tree.root.problem.f(x) - # build the root LMO - free_model(tree.root.problem.tlmo.blmo) - build_LMO( - tree.root.problem.tlmo, - tree.root.problem.integer_variable_bounds, - IntegerBounds(), - tree.root.problem.integer_variables, - ) - gradient = zeros(tree.root.problem.nvars) - tree.root.problem.g(gradient, x) - v = compute_extreme_point(tree.root.problem.tlmo, gradient) - root_dual_gap = FrankWolfe.fast_dot(gradient, x) - FrankWolfe.fast_dot(gradient, v) - - free_model(tree.root.problem.tlmo.blmo) - build_LMO( - tree.root.problem.tlmo, - tree.root.problem.integer_variable_bounds, - node.local_bounds, - tree.root.problem.integer_variables, - ) - - #new_bound_left = M^(- 1 / θ) * (sqrt(new_bound_left) - M * root_dual_gap^θ)^(1 / θ) + fx - root_dual_gap - #new_bound_right = M^(- 1 / θ) * (sqrt(new_bound_right) - M * root_dual_gap^θ)^(1 / θ) + fx - root_dual_gap new_bound_left = M^(- 1 / θ) * (sqrt(new_bound_left) - M * node.dual_gap^θ)^(1 / θ) + fx - node.dual_gap new_bound_right = M^(- 1 / θ) * (sqrt(new_bound_right) - M * node.dual_gap^θ)^(1 / θ) + fx - node.dual_gap From 48ccb45ea0c8b8abc9fc11794c9b754ac5d599b8 Mon Sep 17 00:00:00 2001 From: dhendryc <92737336+dhendryc@users.noreply.github.com> Date: Tue, 3 Sep 2024 14:39:10 +0200 Subject: [PATCH 11/45] Update src/custom_bonobo.jl MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Mathieu Besançon --- src/custom_bonobo.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/custom_bonobo.jl b/src/custom_bonobo.jl index d86094e49..422ba53f8 100644 --- a/src/custom_bonobo.jl +++ b/src/custom_bonobo.jl @@ -75,7 +75,7 @@ function Bonobo.optimize!( Bonobo.branch!(tree, node) callback(tree, node) end - # To make sure that we collect the statistics in case the time limit is crashed. + # To make sure that we collect the statistics in case the time limit is reached. if !haskey(tree.root.result, :global_tightenings) y = Bonobo.get_solution(tree) vertex_storage = FrankWolfe.DeletedVertexStorage(typeof(y)[], 1) From 5650287a9d5e352727c27e149654beb7cf86f289 Mon Sep 17 00:00:00 2001 From: Hendrych Date: Wed, 4 Sep 2024 14:35:07 +0200 Subject: [PATCH 12/45] Add a bit of give to the sharpness constant to avoid cutting of optimal solutions. --- src/tightenings.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/tightenings.jl b/src/tightenings.jl index 2d9ead276..ab6b511b6 100644 --- a/src/tightenings.jl +++ b/src/tightenings.jl @@ -218,7 +218,7 @@ function tightening_lowerbound(tree, node, x, lower_bound) @debug "Using sharpness θ=$θ and M=$M" fx = tree.root.problem.f(x) - sharpness_bound = M^(- 1 / θ) * (sqrt(bound_improvement) - M * node.dual_gap^θ)^(1 / θ) + fx - node.dual_gap + sharpness_bound = M^(- 1 / θ) * 1/ 2 * (sqrt(bound_improvement) - M / 2 * node.dual_gap^θ)^(1 / θ) + fx - node.dual_gap @debug "Sharpness: $lower_bound -> $sharpness_bound" @assert num_fractional == 0 || sharpness_bound > lower_bound end @@ -271,8 +271,8 @@ function prune_children(tree, node, lower_bound_base, x, vidx) elseif M > 0 && θ != Inf fx = tree.root.problem.f(x) - new_bound_left = M^(- 1 / θ) * (sqrt(new_bound_left) - M * node.dual_gap^θ)^(1 / θ) + fx - node.dual_gap - new_bound_right = M^(- 1 / θ) * (sqrt(new_bound_right) - M * node.dual_gap^θ)^(1 / θ) + fx - node.dual_gap + new_bound_left = M^(- 1 / θ) * 1 / 2 * (sqrt(new_bound_left) - M / 2 * node.dual_gap^θ)^(1 / θ) + fx - node.dual_gap + new_bound_right = M^(- 1 / θ) * 1 / 2 * (sqrt(new_bound_right) - M / 2 * node.dual_gap^θ)^(1 / θ) + fx - node.dual_gap if new_bound_left > tree.incumbent @debug "prune left, from $(node.lb) -> $new_bound_left, ub $(tree.incumbent), lb $(node.lb)" From 26b52a08f5f3f04c4a218469c7bba98252ea89f5 Mon Sep 17 00:00:00 2001 From: Hendrych Date: Fri, 6 Sep 2024 10:55:28 +0200 Subject: [PATCH 13/45] Added small definition of sharpness. --- src/interface.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/interface.jl b/src/interface.jl index 920ad881f..7aceb6539 100644 --- a/src/interface.jl +++ b/src/interface.jl @@ -44,6 +44,9 @@ global_dual_tightening - dual tightening maintained globally valid (when new sol 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 sharpness_constant - the constant M > 0 for (θ, M)-sharpness. + f is (θ, M)-sharpness: f satisfies + min_{x* ∈ X*} || x - x* || ≤ M (f(x) - f*) + 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 From 6cd74a73e477d2eb6bc6f769c80c598ce74710f8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Tue, 10 Sep 2024 16:41:23 +0200 Subject: [PATCH 14/45] Update src/callbacks.jl --- src/callbacks.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/callbacks.jl b/src/callbacks.jl index 8da2eff98..ca070317b 100644 --- a/src/callbacks.jl +++ b/src/callbacks.jl @@ -106,8 +106,7 @@ end Branch-and-Bound Callback. Collects statistics and prints logs if verbose is turned on. -Output of Boscia - +Output of Boscia: iter : current iteration of Boscia node id : current node id lower bound : tree_lb(tree) From 7aab85bc0c7c632852752cde1ae659c0c61e6716 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Tue, 10 Sep 2024 16:41:31 +0200 Subject: [PATCH 15/45] Update src/callbacks.jl --- src/callbacks.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/callbacks.jl b/src/callbacks.jl index ca070317b..b803c8ada 100644 --- a/src/callbacks.jl +++ b/src/callbacks.jl @@ -2,7 +2,7 @@ Frank-Wolfe Callback. Is called in every Frank-Wolfe iteration. -Node evaluation can be dymamically stopped here. +Node evaluation can be dynamically stopped here. Time limit is checked. If the vertex is providing a better incumbent, it is added as solution. """ From 32627eef23dac0ba8fc1d7606bdd0dcd272811a0 Mon Sep 17 00:00:00 2001 From: Hendrych Date: Wed, 11 Sep 2024 17:46:35 +0200 Subject: [PATCH 16/45] Don't create child node if the domain oracle is infeasible. --- src/node.jl | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/node.jl b/src/node.jl index c8af87ac3..8daaed9d5 100644 --- a/src/node.jl +++ b/src/node.jl @@ -161,15 +161,18 @@ function Bonobo.get_branching_nodes_info(tree::Bonobo.BnBTree, node::FrankWolfeN x_right = FrankWolfe.compute_active_set_iterate!(active_set_right) domain_oracle = tree.root.options[:domain_oracle] - nodes = if !prune_left && !prune_right + domain_right = domain_oracle(x_right) + domain_left = domain_oracle(x_left) + + nodes = if !prune_left && !prune_right && domain_right && domain_left [node_info_left, node_info_right] elseif prune_left [node_info_right] elseif prune_right [node_info_left] - elseif domain_oracle(x_right) + elseif domain_right #domain_oracle(x_right) [node_info_right] - elseif domain_oracle(x_left) + elseif domain_left #domain_oracle(x_left) [node_info_left] else warn("No childern nodes can be created.") From 414a3f90bf980b741811e05d1ce690578bac52a5 Mon Sep 17 00:00:00 2001 From: Hendrych Date: Wed, 11 Sep 2024 17:53:39 +0200 Subject: [PATCH 17/45] Set version up after bug fix. --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 710426c03..6e1fb7f06 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Boscia" uuid = "36b166db-dac5-4d05-b36a-e6c4cef071c9" authors = ["ZIB IOL"] -version = "0.1.27" +version = "0.1.28" [deps] Bonobo = "f7b14807-3d4d-461a-888a-05dd4bca8bc3" From d79263c814234749b05a2cd46cd22d5455871107 Mon Sep 17 00:00:00 2001 From: Hendrych Date: Thu, 12 Sep 2024 15:53:10 +0200 Subject: [PATCH 18/45] Optimal Experiment Design example. --- examples/oed_utils.jl | 255 ++++++++++++++++++++++++++ examples/optimal_experiment_design.jl | 95 ++++++++++ 2 files changed, 350 insertions(+) create mode 100644 examples/oed_utils.jl create mode 100644 examples/optimal_experiment_design.jl diff --git a/examples/oed_utils.jl b/examples/oed_utils.jl new file mode 100644 index 000000000..575e69b76 --- /dev/null +++ b/examples/oed_utils.jl @@ -0,0 +1,255 @@ +""" + build_data + +seed - for the Random functions. +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(seed, m) + # set up + Random.seed!(seed) + + n = Int(floor(m/10)) + + B = rand(m,n) + B = B'*B + @assert isposdef(B) + D = MvNormal(randn(n),B) + + A = rand(D, m)' + @assert rank(A) == n + + N = floor(1.5*n) + u = floor(N/3) + ub = rand(1.0:u, m) + + return A, n, N, ub +end + +""" +Build Probability Simplex BLMO for Boscia +""" +function build_blmo(m, N, ub) + simplex_lmo = Boscia.ProbabilitySimplexSimpleBLMO(N) + blmo = Boscia.ManagedBoundedLMO(simplex_lmo, fill(0.0, m), ub, collect(1:m), m) + return blmo +end + +""" +Build function and the gradient for the A-criterion. +There is the option to build a safe version of the objective and gradient. If x does not +satisfy the domain oracle, infinity is returned. +If one builds the unsafe version, a FrankWolfe line search must be chosen that can take +a domain oracle as an input like Secant or MonotonicGenericStepSize. +""" +function build_a_criterion(A; μ=1e-4, build_safe=false) + m, n = size(A) + a = m + domain_oracle = build_domain_oracle(A, n) + + function f_a(x) + X = transpose(A)*diagm(x)*A + Matrix(μ *I, n, n) + X = Symmetric(X) + U = cholesky(X) + X_inv = U \ I + return LinearAlgebra.tr(X_inv)/a + end + + function grad_a!(storage, x) + X = transpose(A)*diagm(x)*A + Matrix(μ *I, n, n) + X = Symmetric(X*X) + F = cholesky(X) + for i in 1:length(x) + storage[i] = LinearAlgebra.tr(- (F \ A[i,:]) * transpose(A[i,:]))/a + end + return storage + end + + function f_a_safe(x) + if !domain_oracle(x) + return Inf + end + X = transpose(A)*diagm(x)*A + Matrix(μ *I, n, n) + X = Symmetric(X) + X_inv = LinearAlgebra.inv(X) + return LinearAlgebra.tr(X_inv)/a + end + + function grad_a_safe!(storage, x) + if !domain_oracle(x) + return fill(Inf, length(x)) + end + X = transpose(A)*diagm(x)*A + Matrix(μ *I, n, n) + X = Symmetric(X*X) + F = cholesky(X) + for i in 1:length(x) + storage[i] = LinearAlgebra.tr(- (F \ A[i,:]) * transpose(A[i,:]))/a + end + return storage + end + + if build_safe + return f_a_safe, grad_a_safe! + end + + return f_a, grad_a! +end + +""" +Build function and gradient for the D-criterion. +There is the option to build a safe version of the objective and gradient. If x does not +satisfy the domain oracle, infinity is returned. +If one builds the unsafe version, a FrankWolfe line search must be chosen that can take +a domain oracle as an input like Secant or MonotonicGenericStepSize. +""" +function build_d_criterion(A; μ =0.0, build_safe=false) + m, n = size(A) + a = m + domain_oracle = build_domain_oracle(A, n) + + function f_d(x) + X = transpose(A)*diagm(x)*A + Matrix(μ *I, n, n) + X = Symmetric(X) + return float(-log(det(X))/a) + end + + function grad_d!(storage, x) + X = transpose(A)*diagm(x)*A + Matrix(μ *I, n, n) + X = Symmetric(X) + F = cholesky(X) + for i in 1:length(x) + storage[i] = 1/a * LinearAlgebra.tr(-(F \ A[i,:] )*transpose(A[i,:])) * 1/γ + end + return storage + end + + function f_d_safe(x) + if !domain_oracle(x) + return Inf + end + X = transpose(A)*diagm(x)*A + Matrix(μ *I, n, n) + X = Symmetric(X) + return float(-log(det(X))/a) + end + + function grad_d_safe!(storage, x) + if !domain_oracle(x) + return fill(Inf, length(x)) + end + X = transpose(A)*diagm(x)*A + Matrix(μ *I, n, n) + X = Symmetric(X) + F = cholesky(X) + for i in 1:length(x) + storage[i] = 1/a * LinearAlgebra.tr(-(F \ A[i,:] )*transpose(A[i,:])) * 1/γ + end + # https://stackoverflow.com/questions/46417005/exclude-elements-of-array-based-on-index-julia + return storage + end + + if build_safe + return f_d_safe, grad_d_safe! + end + + return f_d, grad_d! +end + +""" +Find n linearly independent rows of A to build the starting point. +""" +function linearly_independent_rows(A) + S = [] + m, n = size(A) + for i in 1:m + S_i= vcat(S, i) + if rank(A[S_i,:])==length(S_i) + S=S_i + end + if length(S) == n # we only n linearly independent points + return S + end + end + return S # then x= zeros(m) and x[S] = 1 +end + +""" +Build start point used in FrankWolfe and Boscia for the A-Optimal and D-Optimal Design Problem. +The functions are self concordant and so not every point in the feasible region +is in the domain of f and grad!. +""" +function build_start_point(A, N, ub) + # Get n linearly independent rows of A + m, n = size(A) + S = linearly_independent_rows(A,m,n) + @assert length(S) == n + + x = zeros(m) + E = [] + V = Vector{Float64}[] + + while !isempty(setdiff(S, findall(x-> !(iszero(x)),x))) + v = zeros(m) + while sum(v) < N + idx = isempty(setdiff(S, findall(x-> !(iszero(x)),v))) ? rand(setdiff(collect(1:m), S)) : rand(setdiff(S, findall(x-> !(iszero(x)),v))) + if !isapprox(v[idx], 0.0) + @debug "Index $(idx) already picked" + continue + end + v[idx] = min(ub[idx], N - sum(v)) + push!(E, idx) + end + push!(V,v) + x = sum(V .* 1/length(V)) + end + unique!(V) + a = length(V) + x = sum(V .* 1/a) + active_set= FrankWolfe.ActiveSet(fill(1/a, a), V, x) + + return x, active_set, S +end + +""" +Create first incumbent for Boscia and custom BB in a greedy fashion. +""" +function greedy_incumbent(A, N, ub) + # Get n linearly independent rows of A + m, n = size(A) + S = linearly_independent_rows(A,m,n) + @assert length(S) == n + + # set entries to their upper bound + x = zeros(m) + x[S] .= ub[S] + + if isapprox(sum(x), N; atol=1e-4, rtol=1e-2) + return x + elseif sum(x) > N + while sum(x) > N + remove_from_max(x) + end + elseif sum(x) < N + S1 = S + while sum(x) < N + jdx = rand(setdiff(collect(1:m), S1)) + x[jdx] = min(N-sum(x), ub[jdx]) + push!(S1,jdx) + sort!(S1) + end + end + @assert isapprox(sum(x), N; atol=1e-4, rtol=1e-2) + @assert sum(ub - x .>= 0) == m + return x +end + +""" +Check if given point is in the domain of f, i.e. X = transpose(A) * diagm(x) * A +positive definite. +""" +function build_domain_oracle(A, n) + return function domain_oracle(x) + S = findall(x-> !iszero(x),x) + #@show rank(A[S,:]) == n + return size(S) >= n && rank(A[S,:]) == n + end +end \ No newline at end of file diff --git a/examples/optimal_experiment_design.jl b/examples/optimal_experiment_design.jl new file mode 100644 index 000000000..fc18dc38f --- /dev/null +++ b/examples/optimal_experiment_design.jl @@ -0,0 +1,95 @@ +using Boscia +using Random +using Distributions +using LinearAlegbra +using FrankWolfe +using Statistics + +# The function building the problem data and other structures is in a separate file. +include("oed_utils.jl") + +""" + The Optimal Experiment Design Problem consists of choosing a subset of experiments + maximising the information gain. + We generate the Experiment matrix A ∈ R^{mxn} randomly. The information matrix is a linear + map X(x) = A' * diag(x) * A. There exists different information measures Φ. + We concentrate on the A-criterion and D-criterion, i.e. + + Trace(X^{-1}) (A-Criterion) + and + -logdet(X) (D-Criterion). + + Consequently, the optimization problems we want to solve are + + min_x Trace( (A' * diag(x) * A)^{-1} ) + s.t. ∑ x_i = N (A-Optimal Design Problem) + 0 ≤ x ≤ u + + min_x -logdet(A' * diag(x) * A) + s.t. ∑ x_i = N (D-Optimal Design Problem) + 0 ≤ x ≤ u + + where N is our bugdet for the experiments, i.e. this is the amount of experiments + we can perform. We set N = 3/2 * n. The upperbounds u are randomly generated. + + Also, check this paper: https://arxiv.org/abs/2312.11200 and the corresponding + respository https://github.com/ZIB-IOL/OptimalDesignWithBoscia. + + A continuous version of the problem can be found in the examples in FrankWolfe.jl: + https://github.com/ZIB-IOL/FrankWolfe.jl/blob/master/examples/optimal_experiment_design.jl +""" + +seed = 1234 +m = 50 + +## A-Optimal Design Problem + +A, N, ub = build_data(seed, m) + +# sharpness constants +σ = minimum(A' * A) +a = maximum(ub) * maximum([norm(A[i,:])^2 for i=1:size(A,1)]) +θ = 1/2 +M = sqrt(2) * a / (σ^2 * sqrt(2 * n)) + +f, grad! = build_a_criterion(A, build_safe=true) +blmo = build_blmo(m, N, ub) +x0, active_set = build_start_point(A, N, ub) +z = build_greedy_incumbent(A, N, ub) + +x, _ result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, sharpness_exponent=θ, sharpness_constant=M) + +f, grad! = build_a_criterion(A, build_safe=false) +blmo = build_blmo(m, N, ub) +x0, active_set = build_start_point(A, N, ub) +z = build_greedy_incumbent(A, N, ub) +domain_oracle = build_domain_oracle(A, n) + +x, _ result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, line_search=FrankWolfe.Secant(40, 1e-8, domain_oracle), sharpness_exponent=θ, sharpness_constant=M) + + +## D-Optimal Design Problem + +A, N, ub = build_data(seed, m) + +# sharpness constants +σ = minimum(A' * A) +a = maximum(ub) * maximum([norm(A[i,:])^2 for i=1:size(A,1)]) +θ = 1/2 +M = sqrt(2) * a / (σ^2 * sqrt(n)) + +f, grad! = build_d_criterion(A, build_safe=true) +blmo = build_blmo(m, N, ub) +x0, active_set = build_start_point(A, N, ub) +z = build_greedy_incumbent(A, N, ub) + +x, _ result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, sharpness_exponent=θ, sharpness_constant=M) + +f, grad! = build_d_criterion(A, build_safe=false) +blmo = build_blmo(m, N, ub) +x0, active_set = build_start_point(A, N, ub) +z = build_greedy_incumbent(A, N, ub) +domain_oracle = build_domain_oracle(A, n) + +x, _ result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, line_search=FrankWolfe.Secant(40, 1e-8, domain_oracle), sharpness_exponent=θ, sharpness_constant=M) + From 802f3b8cc6bee27b78b78be45ae998a842e536d4 Mon Sep 17 00:00:00 2001 From: Hendrych Date: Thu, 12 Sep 2024 16:00:13 +0200 Subject: [PATCH 19/45] Fix syntax issues. --- examples/oed_utils.jl | 36 ++++++++++++++++++++++++--- examples/optimal_experiment_design.jl | 22 ++++++++-------- 2 files changed, 44 insertions(+), 14 deletions(-) diff --git a/examples/oed_utils.jl b/examples/oed_utils.jl index 575e69b76..42ae65837 100644 --- a/examples/oed_utils.jl +++ b/examples/oed_utils.jl @@ -172,6 +172,36 @@ function linearly_independent_rows(A) return S # then x= zeros(m) and x[S] = 1 end +function add_to_min(x, u) + perm = sortperm(x) + j = findfirst(x->x != 0, x[perm]) + + for i in j:length(x) + if x[perm[i]] < u[perm[i]] + x[perm[i]] += 1 + break + else + continue + end + end + return x +end + +function remove_from_max(x) + perm = sortperm(x, rev = true) + j = findlast(x->x != 0, x[perm]) + + for i in 1:j + if x[perm[i]] > 1 + x[perm[i]] -= 1 + break + else + continue + end + end + return x +end + """ Build start point used in FrankWolfe and Boscia for the A-Optimal and D-Optimal Design Problem. The functions are self concordant and so not every point in the feasible region @@ -180,7 +210,7 @@ is in the domain of f and grad!. function build_start_point(A, N, ub) # Get n linearly independent rows of A m, n = size(A) - S = linearly_independent_rows(A,m,n) + S = linearly_independent_rows(A) @assert length(S) == n x = zeros(m) @@ -215,7 +245,7 @@ Create first incumbent for Boscia and custom BB in a greedy fashion. function greedy_incumbent(A, N, ub) # Get n linearly independent rows of A m, n = size(A) - S = linearly_independent_rows(A,m,n) + S = linearly_independent_rows(A) @assert length(S) == n # set entries to their upper bound @@ -250,6 +280,6 @@ function build_domain_oracle(A, n) return function domain_oracle(x) S = findall(x-> !iszero(x),x) #@show rank(A[S,:]) == n - return size(S) >= n && rank(A[S,:]) == n + return length(S) >= n && rank(A[S,:]) == n end end \ No newline at end of file diff --git a/examples/optimal_experiment_design.jl b/examples/optimal_experiment_design.jl index fc18dc38f..cc3641ac1 100644 --- a/examples/optimal_experiment_design.jl +++ b/examples/optimal_experiment_design.jl @@ -1,7 +1,7 @@ using Boscia using Random using Distributions -using LinearAlegbra +using LinearAlgebra using FrankWolfe using Statistics @@ -44,7 +44,7 @@ m = 50 ## A-Optimal Design Problem -A, N, ub = build_data(seed, m) +A, n, N, ub = build_data(seed, m) # sharpness constants σ = minimum(A' * A) @@ -55,22 +55,22 @@ M = sqrt(2) * a / (σ^2 * sqrt(2 * n)) f, grad! = build_a_criterion(A, build_safe=true) blmo = build_blmo(m, N, ub) x0, active_set = build_start_point(A, N, ub) -z = build_greedy_incumbent(A, N, ub) +z = greedy_incumbent(A, N, ub) -x, _ result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, sharpness_exponent=θ, sharpness_constant=M) +x, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=true, sharpness_exponent=θ, sharpness_constant=M) f, grad! = build_a_criterion(A, build_safe=false) blmo = build_blmo(m, N, ub) x0, active_set = build_start_point(A, N, ub) -z = build_greedy_incumbent(A, N, ub) +z = greedy_incumbent(A, N, ub) domain_oracle = build_domain_oracle(A, n) -x, _ result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, line_search=FrankWolfe.Secant(40, 1e-8, domain_oracle), sharpness_exponent=θ, sharpness_constant=M) +x, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=true, line_search=FrankWolfe.Secant(40, 1e-8, domain_oracle), sharpness_exponent=θ, sharpness_constant=M) ## D-Optimal Design Problem -A, N, ub = build_data(seed, m) +A, n, N, ub = build_data(seed, m) # sharpness constants σ = minimum(A' * A) @@ -81,15 +81,15 @@ M = sqrt(2) * a / (σ^2 * sqrt(n)) f, grad! = build_d_criterion(A, build_safe=true) blmo = build_blmo(m, N, ub) x0, active_set = build_start_point(A, N, ub) -z = build_greedy_incumbent(A, N, ub) +z = greedy_incumbent(A, N, ub) -x, _ result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, sharpness_exponent=θ, sharpness_constant=M) +x, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=true, sharpness_exponent=θ, sharpness_constant=M) f, grad! = build_d_criterion(A, build_safe=false) blmo = build_blmo(m, N, ub) x0, active_set = build_start_point(A, N, ub) -z = build_greedy_incumbent(A, N, ub) +z = greedy_incumbent(A, N, ub) domain_oracle = build_domain_oracle(A, n) -x, _ result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, line_search=FrankWolfe.Secant(40, 1e-8, domain_oracle), sharpness_exponent=θ, sharpness_constant=M) +x, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=true, line_search=FrankWolfe.Secant(40, 1e-8, domain_oracle), sharpness_exponent=θ, sharpness_constant=M) From c99356238cbac4000e61d6ff1e98f5e5ca0ec22a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Fri, 13 Sep 2024 09:33:03 +0200 Subject: [PATCH 20/45] Apply suggestions from code review --- src/node.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/node.jl b/src/node.jl index 8daaed9d5..5bb18bcea 100644 --- a/src/node.jl +++ b/src/node.jl @@ -170,9 +170,9 @@ function Bonobo.get_branching_nodes_info(tree::Bonobo.BnBTree, node::FrankWolfeN [node_info_right] elseif prune_right [node_info_left] - elseif domain_right #domain_oracle(x_right) + elseif domain_right # x_right in domain [node_info_right] - elseif domain_left #domain_oracle(x_left) + elseif domain_left # x_left in domain [node_info_left] else warn("No childern nodes can be created.") From a068fc8c4407c50c77d7467a72f3a44700bf7725 Mon Sep 17 00:00:00 2001 From: Hendrych Date: Fri, 13 Sep 2024 12:56:13 +0200 Subject: [PATCH 21/45] Add corrected sharpness constants. --- examples/oed_utils.jl | 4 ++-- examples/optimal_experiment_design.jl | 15 +++++++++------ 2 files changed, 11 insertions(+), 8 deletions(-) diff --git a/examples/oed_utils.jl b/examples/oed_utils.jl index 42ae65837..d81580fb2 100644 --- a/examples/oed_utils.jl +++ b/examples/oed_utils.jl @@ -119,7 +119,7 @@ function build_d_criterion(A; μ =0.0, build_safe=false) X = Symmetric(X) F = cholesky(X) for i in 1:length(x) - storage[i] = 1/a * LinearAlgebra.tr(-(F \ A[i,:] )*transpose(A[i,:])) * 1/γ + storage[i] = 1/a * LinearAlgebra.tr(-(F \ A[i,:] )*transpose(A[i,:])) end return storage end @@ -141,7 +141,7 @@ function build_d_criterion(A; μ =0.0, build_safe=false) X = Symmetric(X) F = cholesky(X) for i in 1:length(x) - storage[i] = 1/a * LinearAlgebra.tr(-(F \ A[i,:] )*transpose(A[i,:])) * 1/γ + storage[i] = 1/a * LinearAlgebra.tr(-(F \ A[i,:] )*transpose(A[i,:])) end # https://stackoverflow.com/questions/46417005/exclude-elements-of-array-based-on-index-julia return storage diff --git a/examples/optimal_experiment_design.jl b/examples/optimal_experiment_design.jl index cc3641ac1..5dd4ca379 100644 --- a/examples/optimal_experiment_design.jl +++ b/examples/optimal_experiment_design.jl @@ -41,6 +41,7 @@ include("oed_utils.jl") seed = 1234 m = 50 +verbose = true ## A-Optimal Design Problem @@ -50,14 +51,15 @@ A, n, N, ub = build_data(seed, m) σ = minimum(A' * A) a = maximum(ub) * maximum([norm(A[i,:])^2 for i=1:size(A,1)]) θ = 1/2 -M = sqrt(2) * a / (σ^2 * sqrt(2 * n)) +M = n * σ^4 / a^2 + f, grad! = build_a_criterion(A, build_safe=true) blmo = build_blmo(m, N, ub) x0, active_set = build_start_point(A, N, ub) z = greedy_incumbent(A, N, ub) -x, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=true, sharpness_exponent=θ, sharpness_constant=M) +x, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, sharpness_exponent=θ, sharpness_constant=M) f, grad! = build_a_criterion(A, build_safe=false) blmo = build_blmo(m, N, ub) @@ -65,7 +67,7 @@ x0, active_set = build_start_point(A, N, ub) z = greedy_incumbent(A, N, ub) domain_oracle = build_domain_oracle(A, n) -x, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=true, line_search=FrankWolfe.Secant(40, 1e-8, domain_oracle), sharpness_exponent=θ, sharpness_constant=M) +x, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, line_search=FrankWolfe.Secant(40, 1e-8, domain_oracle), sharpness_exponent=θ, sharpness_constant=M) ## D-Optimal Design Problem @@ -76,14 +78,15 @@ A, n, N, ub = build_data(seed, m) σ = minimum(A' * A) a = maximum(ub) * maximum([norm(A[i,:])^2 for i=1:size(A,1)]) θ = 1/2 -M = sqrt(2) * a / (σ^2 * sqrt(n)) +M = n * σ^4 / (2 * a^2) + f, grad! = build_d_criterion(A, build_safe=true) blmo = build_blmo(m, N, ub) x0, active_set = build_start_point(A, N, ub) z = greedy_incumbent(A, N, ub) -x, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=true, sharpness_exponent=θ, sharpness_constant=M) +x, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, sharpness_exponent=θ, sharpness_constant=M) f, grad! = build_d_criterion(A, build_safe=false) blmo = build_blmo(m, N, ub) @@ -91,5 +94,5 @@ x0, active_set = build_start_point(A, N, ub) z = greedy_incumbent(A, N, ub) domain_oracle = build_domain_oracle(A, n) -x, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=true, line_search=FrankWolfe.Secant(40, 1e-8, domain_oracle), sharpness_exponent=θ, sharpness_constant=M) +x, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, line_search=FrankWolfe.Secant(40, 1e-8, domain_oracle), sharpness_exponent=θ, sharpness_constant=M, fw_verbose=true) From b1b52736722abcc28157ae43fc8697a7a16ada32 Mon Sep 17 00:00:00 2001 From: Hendrych Date: Fri, 13 Sep 2024 15:22:47 +0200 Subject: [PATCH 22/45] In case the vertex is domain infeasible, use x to get the data type of f(x). --- src/interface.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/interface.jl b/src/interface.jl index 5a49b70bd..36498de5b 100644 --- a/src/interface.jl +++ b/src/interface.jl @@ -159,8 +159,8 @@ function solve( for a in active_set.atoms @assert is_linear_feasible(blmo, a) end - v = active_set.atoms[1] x = FrankWolfe.compute_active_set_iterate!(active_set) + v = x @assert isfinite(f(x)) end vertex_storage = FrankWolfe.DeletedVertexStorage(typeof(v)[], 1) From b676ef6452f41527d8382c936938a37717af9cbd Mon Sep 17 00:00:00 2001 From: Hendrych Date: Fri, 13 Sep 2024 15:37:13 +0200 Subject: [PATCH 23/45] Syntax issue fix. --- src/node.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/node.jl b/src/node.jl index 3d1915cd8..d6b5b9323 100644 --- a/src/node.jl +++ b/src/node.jl @@ -175,7 +175,7 @@ function Bonobo.get_branching_nodes_info(tree::Bonobo.BnBTree, node::FrankWolfeN elseif domain_left #domain_oracle(x_left) [node_info_left] else - warn("No childern nodes can be created.") + @warn "No childern nodes can be created." end return nodes end From 1f26759e98cacf119c72f60df85b7d4d058645aa Mon Sep 17 00:00:00 2001 From: Hendrych Date: Fri, 13 Sep 2024 15:37:28 +0200 Subject: [PATCH 24/45] Corrections. --- examples/oed_utils.jl | 16 +++++++++++++--- examples/optimal_experiment_design.jl | 15 +++++++++------ 2 files changed, 22 insertions(+), 9 deletions(-) diff --git a/examples/oed_utils.jl b/examples/oed_utils.jl index d81580fb2..b67bfb42b 100644 --- a/examples/oed_utils.jl +++ b/examples/oed_utils.jl @@ -105,7 +105,7 @@ a domain oracle as an input like Secant or MonotonicGenericStepSize. """ function build_d_criterion(A; μ =0.0, build_safe=false) m, n = size(A) - a = m + a = 1#m domain_oracle = build_domain_oracle(A, n) function f_d(x) @@ -275,11 +275,21 @@ end """ Check if given point is in the domain of f, i.e. X = transpose(A) * diagm(x) * A positive definite. + +(a) Check the rank of A restricted to the rows activated by x. +(b) Check if the resulting information matrix A' * diagm(x) * A is psd. + +(b) is a bit faster for smaller dimensions (< 100). For larger (> 200) (a) is faster. """ function build_domain_oracle(A, n) return function domain_oracle(x) S = findall(x-> !iszero(x),x) - #@show rank(A[S,:]) == n return length(S) >= n && rank(A[S,:]) == n end -end \ No newline at end of file +end + +function build_domain_oracle2(A) + return function domain_oracle2(x) + return isposdef(Symmetric(A' * diagm(x) * A)) + end +end diff --git a/examples/optimal_experiment_design.jl b/examples/optimal_experiment_design.jl index 5dd4ca379..64841d37d 100644 --- a/examples/optimal_experiment_design.jl +++ b/examples/optimal_experiment_design.jl @@ -40,7 +40,7 @@ include("oed_utils.jl") """ seed = 1234 -m = 50 +m = 100 verbose = true ## A-Optimal Design Problem @@ -58,8 +58,9 @@ f, grad! = build_a_criterion(A, build_safe=true) blmo = build_blmo(m, N, ub) x0, active_set = build_start_point(A, N, ub) z = greedy_incumbent(A, N, ub) +domain_oracle = build_domain_oracle(A, n) -x, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, sharpness_exponent=θ, sharpness_constant=M) +x, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, sharpness_exponent=θ, sharpness_constant=M, domain_oracle=domain_oracle) f, grad! = build_a_criterion(A, build_safe=false) blmo = build_blmo(m, N, ub) @@ -67,7 +68,7 @@ x0, active_set = build_start_point(A, N, ub) z = greedy_incumbent(A, N, ub) domain_oracle = build_domain_oracle(A, n) -x, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, line_search=FrankWolfe.Secant(40, 1e-8, domain_oracle), sharpness_exponent=θ, sharpness_constant=M) +x, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, line_search=FrankWolfe.Secant(40, 1e-8, domain_oracle), sharpness_exponent=θ, sharpness_constant=M, domain_oracle=domain_oracle) ## D-Optimal Design Problem @@ -85,8 +86,11 @@ f, grad! = build_d_criterion(A, build_safe=true) blmo = build_blmo(m, N, ub) x0, active_set = build_start_point(A, N, ub) z = greedy_incumbent(A, N, ub) +domain_oracle = build_domain_oracle(A, n) -x, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, sharpness_exponent=θ, sharpness_constant=M) +x, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, sharpness_exponent=θ, sharpness_constant=M, domain_oracle=domain_oracle) + +A, n, N, ub = build_data(seed, m) f, grad! = build_d_criterion(A, build_safe=false) blmo = build_blmo(m, N, ub) @@ -94,5 +98,4 @@ x0, active_set = build_start_point(A, N, ub) z = greedy_incumbent(A, N, ub) domain_oracle = build_domain_oracle(A, n) -x, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, line_search=FrankWolfe.Secant(40, 1e-8, domain_oracle), sharpness_exponent=θ, sharpness_constant=M, fw_verbose=true) - +x, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, line_search=FrankWolfe.Secant(40, 1e-8, domain_oracle), sharpness_exponent=θ, sharpness_constant=M, domain_oracle=domain_oracle) From 98783d53141d34f239aba73f028dd864265a7960 Mon Sep 17 00:00:00 2001 From: Hendrych Date: Fri, 13 Sep 2024 15:51:39 +0200 Subject: [PATCH 25/45] Print statement in case assert fails. --- src/interface.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/interface.jl b/src/interface.jl index b496a1ea8..74d0f30fe 100644 --- a/src/interface.jl +++ b/src/interface.jl @@ -642,7 +642,8 @@ function postsolve(tree, result, time_ref, verbose, max_iteration_post) else @info "primal >= tree.incumbent" @assert primal <= tree.incumbent + 1e-3 || - isapprox(primal, tree.incumbent, atol=1e-6, rtol=1e-2) + isapprox(primal, tree.incumbent, atol=1e-6, rtol=1e-2) "primal <= tree.incumbent + 1e-3 || + isapprox(primal, tree.incumbent, atol=1e-6, rtol=1e-2): primal=$(primal) and tree.incumbent=$(tree.incumbent)" end @info "postsolve did not improve the solution" primal = tree.incumbent_solution.objective = tree.solutions[1].objective From c55223f554d0aa61891dbdbd1627205a7e3e4241 Mon Sep 17 00:00:00 2001 From: Hendrych Date: Mon, 16 Sep 2024 12:11:30 +0200 Subject: [PATCH 26/45] Print logs for better bug hunting. --- test/heuristics.jl | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/test/heuristics.jl b/test/heuristics.jl index 6005832c3..e06d74ffa 100644 --- a/test/heuristics.jl +++ b/test/heuristics.jl @@ -11,6 +11,11 @@ using Dates const MOI = MathOptInterface const MOIU = MOI.Utilities +seed = rand(UInt64) +seed = 0x61746adc8587896d +@show seed +Random.seed!(seed) + n = 20 x_sol = rand(1:floor(Int, n/4), n) N = sum(x_sol) @@ -159,7 +164,7 @@ diffi = x_sol + 0.3 * dir heu = Boscia.Heuristic(Boscia.probability_rounding, 0.6, :probability_rounding) x, _, result = - Boscia.solve(f, grad!, sblmo, fill(0.0, m), fill(1.0, m), int_vars, n, custom_heuristics=[heu], rounding_prob=0.0) + Boscia.solve(f, grad!, sblmo, fill(0.0, m), fill(1.0, m), int_vars, n, custom_heuristics=[heu], rounding_prob=0.0, verbose=true) @test f(x) ≥ f(x_sol) if isapprox(sum(x_sol), N) From b5e7e240f9a239067604476857a369930ee65373 Mon Sep 17 00:00:00 2001 From: Hendrych Date: Mon, 16 Sep 2024 12:54:50 +0200 Subject: [PATCH 27/45] Another try to reproduce error. --- test/heuristics.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/heuristics.jl b/test/heuristics.jl index e06d74ffa..e3daac936 100644 --- a/test/heuristics.jl +++ b/test/heuristics.jl @@ -164,7 +164,7 @@ diffi = x_sol + 0.3 * dir heu = Boscia.Heuristic(Boscia.probability_rounding, 0.6, :probability_rounding) x, _, result = - Boscia.solve(f, grad!, sblmo, fill(0.0, m), fill(1.0, m), int_vars, n, custom_heuristics=[heu], rounding_prob=0.0, verbose=true) + Boscia.solve(f, grad!, sblmo, fill(0.0, m), fill(1.0, m), int_vars, n, custom_heuristics=[heu], rounding_prob=0.0, verbose=false) @test f(x) ≥ f(x_sol) if isapprox(sum(x_sol), N) From 722cb2d65743d8a13cb87b232fe8bcf789aeeabb Mon Sep 17 00:00:00 2001 From: Hendrych Date: Mon, 16 Sep 2024 13:55:47 +0200 Subject: [PATCH 28/45] Test lower bound against tree.incumbent plus current fw_epsilon. --- src/node.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/node.jl b/src/node.jl index 5bb18bcea..170960ebd 100644 --- a/src/node.jl +++ b/src/node.jl @@ -323,7 +323,7 @@ function Bonobo.evaluate_node!(tree::Bonobo.BnBTree, node::FrankWolfeNode) elseif node.id == 1 @debug "Lower bound of root node: $(lower_bound)" @debug "Current incumbent: $(tree.incumbent)" - @assert lower_bound <= tree.incumbent + 1e-5 "lower_bound <= tree.incumbent + 1e-5 : $(lower_bound) <= $(tree.incumbent + 1e-5)" + @assert lower_bound <= tree.incumbent + node.fw_dual_gap_limit "lower_bound <= tree.incumbent + node.fw_dual_gap_limit : $(lower_bound) <= $(tree.incumbent + node.fw_dual_gap_limit)" end # Call heuristic From c6aaf84a9e2a306a04fadaaf59b45154c7401163 Mon Sep 17 00:00:00 2001 From: Hendrych Date: Mon, 16 Sep 2024 14:25:44 +0200 Subject: [PATCH 29/45] Have the same seed as in runtests.jl. --- test/heuristics.jl | 5 ----- 1 file changed, 5 deletions(-) diff --git a/test/heuristics.jl b/test/heuristics.jl index e3daac936..e2422c402 100644 --- a/test/heuristics.jl +++ b/test/heuristics.jl @@ -11,11 +11,6 @@ using Dates const MOI = MathOptInterface const MOIU = MOI.Utilities -seed = rand(UInt64) -seed = 0x61746adc8587896d -@show seed -Random.seed!(seed) - n = 20 x_sol = rand(1:floor(Int, n/4), n) N = sum(x_sol) From dc659622658252480bed0126b2d8b053d26b6316 Mon Sep 17 00:00:00 2001 From: Hendrych Date: Mon, 16 Sep 2024 15:12:30 +0200 Subject: [PATCH 30/45] Rename A to Ex_mat to avoid constant redefinition. --- examples/optimal_experiment_design.jl | 46 +++++++++++++-------------- 1 file changed, 23 insertions(+), 23 deletions(-) diff --git a/examples/optimal_experiment_design.jl b/examples/optimal_experiment_design.jl index 64841d37d..28e5cc0df 100644 --- a/examples/optimal_experiment_design.jl +++ b/examples/optimal_experiment_design.jl @@ -45,57 +45,57 @@ verbose = true ## A-Optimal Design Problem -A, n, N, ub = build_data(seed, m) +Ex_mat, n, N, ub = build_data(seed, m) # sharpness constants -σ = minimum(A' * A) -a = maximum(ub) * maximum([norm(A[i,:])^2 for i=1:size(A,1)]) +σ = minimum(Ex_mat' * Ex_mat) +a = maximum(ub) * maximum([norm(Ex_mat[i,:])^2 for i=1:size(Ex_mat,1)]) θ = 1/2 M = n * σ^4 / a^2 -f, grad! = build_a_criterion(A, build_safe=true) +f, grad! = build_a_criterion(Ex_mat, build_safe=true) blmo = build_blmo(m, N, ub) -x0, active_set = build_start_point(A, N, ub) -z = greedy_incumbent(A, N, ub) -domain_oracle = build_domain_oracle(A, n) +x0, active_set = build_start_point(Ex_mat, N, ub) +z = greedy_incumbent(Ex_mat, N, ub) +domain_oracle = build_domain_oracle(Ex_mat, n) x, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, sharpness_exponent=θ, sharpness_constant=M, domain_oracle=domain_oracle) -f, grad! = build_a_criterion(A, build_safe=false) +f, grad! = build_a_criterion(Ex_mat, build_safe=false) blmo = build_blmo(m, N, ub) -x0, active_set = build_start_point(A, N, ub) -z = greedy_incumbent(A, N, ub) -domain_oracle = build_domain_oracle(A, n) +x0, active_set = build_start_point(Ex_mat, N, ub) +z = greedy_incumbent(Ex_mat, N, ub) +domain_oracle = build_domain_oracle(Ex_mat, n) x, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, line_search=FrankWolfe.Secant(40, 1e-8, domain_oracle), sharpness_exponent=θ, sharpness_constant=M, domain_oracle=domain_oracle) ## D-Optimal Design Problem -A, n, N, ub = build_data(seed, m) +Ex_mat, n, N, ub = build_data(seed, m) # sharpness constants -σ = minimum(A' * A) -a = maximum(ub) * maximum([norm(A[i,:])^2 for i=1:size(A,1)]) +σ = minimum(Ex_mat' * Ex_mat) +a = maximum(ub) * maximum([norm(Ex_mat[i,:])^2 for i=1:size(Ex_mat,1)]) θ = 1/2 M = n * σ^4 / (2 * a^2) -f, grad! = build_d_criterion(A, build_safe=true) +f, grad! = build_d_criterion(Ex_mat, build_safe=true) blmo = build_blmo(m, N, ub) -x0, active_set = build_start_point(A, N, ub) -z = greedy_incumbent(A, N, ub) -domain_oracle = build_domain_oracle(A, n) +x0, active_set = build_start_point(Ex_mat, N, ub) +z = greedy_incumbent(Ex_mat, N, ub) +domain_oracle = build_domain_oracle(Ex_mat, n) x, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, sharpness_exponent=θ, sharpness_constant=M, domain_oracle=domain_oracle) -A, n, N, ub = build_data(seed, m) +Ex_mat, n, N, ub = build_data(seed, m) -f, grad! = build_d_criterion(A, build_safe=false) +f, grad! = build_d_criterion(Ex_mat, build_safe=false) blmo = build_blmo(m, N, ub) -x0, active_set = build_start_point(A, N, ub) -z = greedy_incumbent(A, N, ub) -domain_oracle = build_domain_oracle(A, n) +x0, active_set = build_start_point(Ex_mat, N, ub) +z = greedy_incumbent(Ex_mat, N, ub) +domain_oracle = build_domain_oracle(Ex_mat, n) x, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, line_search=FrankWolfe.Secant(40, 1e-8, domain_oracle), sharpness_exponent=θ, sharpness_constant=M, domain_oracle=domain_oracle) From 4c16b12f881c325b4d1771c548f2fdf678e5e04c Mon Sep 17 00:00:00 2001 From: Hendrych Date: Mon, 16 Sep 2024 17:39:51 +0200 Subject: [PATCH 31/45] Clean up another constant redefinition. --- examples/optimal_experiment_design.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/optimal_experiment_design.jl b/examples/optimal_experiment_design.jl index 28e5cc0df..1d43ff630 100644 --- a/examples/optimal_experiment_design.jl +++ b/examples/optimal_experiment_design.jl @@ -49,9 +49,9 @@ Ex_mat, n, N, ub = build_data(seed, m) # sharpness constants σ = minimum(Ex_mat' * Ex_mat) -a = maximum(ub) * maximum([norm(Ex_mat[i,:])^2 for i=1:size(Ex_mat,1)]) +λ_max = maximum(ub) * maximum([norm(Ex_mat[i,:])^2 for i=1:size(Ex_mat,1)]) θ = 1/2 -M = n * σ^4 / a^2 +M = n * σ^4 / λ_max^2 f, grad! = build_a_criterion(Ex_mat, build_safe=true) @@ -77,9 +77,9 @@ Ex_mat, n, N, ub = build_data(seed, m) # sharpness constants σ = minimum(Ex_mat' * Ex_mat) -a = maximum(ub) * maximum([norm(Ex_mat[i,:])^2 for i=1:size(Ex_mat,1)]) +λ_max = maximum(ub) * maximum([norm(Ex_mat[i,:])^2 for i=1:size(Ex_mat,1)]) θ = 1/2 -M = n * σ^4 / (2 * a^2) +M = n * σ^4 / (2 * λ_max^2) f, grad! = build_d_criterion(Ex_mat, build_safe=true) From 3a4da67c049b9d592e0555f79094fee9f508a033 Mon Sep 17 00:00:00 2001 From: Hendrych Date: Mon, 16 Sep 2024 18:10:17 +0200 Subject: [PATCH 32/45] Add test wrapper. --- examples/optimal_experiment_design.jl | 83 +++++++++++++++------------ 1 file changed, 46 insertions(+), 37 deletions(-) diff --git a/examples/optimal_experiment_design.jl b/examples/optimal_experiment_design.jl index 1d43ff630..341eefe32 100644 --- a/examples/optimal_experiment_design.jl +++ b/examples/optimal_experiment_design.jl @@ -4,6 +4,7 @@ using Distributions using LinearAlgebra using FrankWolfe using Statistics +using Test # The function building the problem data and other structures is in a separate file. include("oed_utils.jl") @@ -45,57 +46,65 @@ verbose = true ## A-Optimal Design Problem -Ex_mat, n, N, ub = build_data(seed, m) +@testset "A-Optimal Design" begin -# sharpness constants -σ = minimum(Ex_mat' * Ex_mat) -λ_max = maximum(ub) * maximum([norm(Ex_mat[i,:])^2 for i=1:size(Ex_mat,1)]) -θ = 1/2 -M = n * σ^4 / λ_max^2 + Ex_mat, n, N, ub = build_data(seed, m) + # sharpness constants + σ = minimum(Ex_mat' * Ex_mat) + λ_max = maximum(ub) * maximum([norm(Ex_mat[i,:])^2 for i=1:size(Ex_mat,1)]) + θ = 1/2 + M = n * σ^4 / λ_max^2 -f, grad! = build_a_criterion(Ex_mat, build_safe=true) -blmo = build_blmo(m, N, ub) -x0, active_set = build_start_point(Ex_mat, N, ub) -z = greedy_incumbent(Ex_mat, N, ub) -domain_oracle = build_domain_oracle(Ex_mat, n) -x, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, sharpness_exponent=θ, sharpness_constant=M, domain_oracle=domain_oracle) + f, grad! = build_a_criterion(Ex_mat, build_safe=true) + blmo = build_blmo(m, N, ub) + x0, active_set = build_start_point(Ex_mat, N, ub) + z = greedy_incumbent(Ex_mat, N, ub) + domain_oracle = build_domain_oracle(Ex_mat, n) -f, grad! = build_a_criterion(Ex_mat, build_safe=false) -blmo = build_blmo(m, N, ub) -x0, active_set = build_start_point(Ex_mat, N, ub) -z = greedy_incumbent(Ex_mat, N, ub) -domain_oracle = build_domain_oracle(Ex_mat, n) + x, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, sharpness_exponent=θ, sharpness_constant=M, domain_oracle=domain_oracle) -x, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, line_search=FrankWolfe.Secant(40, 1e-8, domain_oracle), sharpness_exponent=θ, sharpness_constant=M, domain_oracle=domain_oracle) + f, grad! = build_a_criterion(Ex_mat, build_safe=false) + blmo = build_blmo(m, N, ub) + x0, active_set = build_start_point(Ex_mat, N, ub) + z = greedy_incumbent(Ex_mat, N, ub) + domain_oracle = build_domain_oracle(Ex_mat, n) + x_s, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, line_search=FrankWolfe.Secant(40, 1e-8, domain_oracle), sharpness_exponent=θ, sharpness_constant=M, domain_oracle=domain_oracle) + + @test isapprox(f(x), f(x_s)) +end ## D-Optimal Design Problem -Ex_mat, n, N, ub = build_data(seed, m) +@testset "D-optimal Design" begin + Ex_mat, n, N, ub = build_data(seed, m) + + # sharpness constants + σ = minimum(Ex_mat' * Ex_mat) + λ_max = maximum(ub) * maximum([norm(Ex_mat[i,:])^2 for i=1:size(Ex_mat,1)]) + θ = 1/2 + M = n * σ^4 / (2 * λ_max^2) -# sharpness constants -σ = minimum(Ex_mat' * Ex_mat) -λ_max = maximum(ub) * maximum([norm(Ex_mat[i,:])^2 for i=1:size(Ex_mat,1)]) -θ = 1/2 -M = n * σ^4 / (2 * λ_max^2) + f, grad! = build_d_criterion(Ex_mat, build_safe=true) + blmo = build_blmo(m, N, ub) + x0, active_set = build_start_point(Ex_mat, N, ub) + z = greedy_incumbent(Ex_mat, N, ub) + domain_oracle = build_domain_oracle(Ex_mat, n) -f, grad! = build_d_criterion(Ex_mat, build_safe=true) -blmo = build_blmo(m, N, ub) -x0, active_set = build_start_point(Ex_mat, N, ub) -z = greedy_incumbent(Ex_mat, N, ub) -domain_oracle = build_domain_oracle(Ex_mat, n) + x, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, sharpness_exponent=θ, sharpness_constant=M, domain_oracle=domain_oracle) -x, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, sharpness_exponent=θ, sharpness_constant=M, domain_oracle=domain_oracle) + Ex_mat, n, N, ub = build_data(seed, m) -Ex_mat, n, N, ub = build_data(seed, m) + f, grad! = build_d_criterion(Ex_mat, build_safe=false) + blmo = build_blmo(m, N, ub) + x0, active_set = build_start_point(Ex_mat, N, ub) + z = greedy_incumbent(Ex_mat, N, ub) + domain_oracle = build_domain_oracle(Ex_mat, n) -f, grad! = build_d_criterion(Ex_mat, build_safe=false) -blmo = build_blmo(m, N, ub) -x0, active_set = build_start_point(Ex_mat, N, ub) -z = greedy_incumbent(Ex_mat, N, ub) -domain_oracle = build_domain_oracle(Ex_mat, n) + x_s, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, line_search=FrankWolfe.Secant(40, 1e-8, domain_oracle), sharpness_exponent=θ, sharpness_constant=M, domain_oracle=domain_oracle) -x, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, line_search=FrankWolfe.Secant(40, 1e-8, domain_oracle), sharpness_exponent=θ, sharpness_constant=M, domain_oracle=domain_oracle) + @test isapprox(f(x), f(x_s)) +end From 015b4a75c6e7ddfb08a7bba4d77138bde2d8693a Mon Sep 17 00:00:00 2001 From: Hendrych Date: Tue, 17 Sep 2024 11:20:13 +0200 Subject: [PATCH 33/45] Debug statements. --- examples/optimal_experiment_design.jl | 26 +++++++++++++++++--------- 1 file changed, 17 insertions(+), 9 deletions(-) diff --git a/examples/optimal_experiment_design.jl b/examples/optimal_experiment_design.jl index 341eefe32..5a36e38bd 100644 --- a/examples/optimal_experiment_design.jl +++ b/examples/optimal_experiment_design.jl @@ -41,7 +41,7 @@ include("oed_utils.jl") """ seed = 1234 -m = 100 +m = 40 verbose = true ## A-Optimal Design Problem @@ -62,18 +62,21 @@ verbose = true x0, active_set = build_start_point(Ex_mat, N, ub) z = greedy_incumbent(Ex_mat, N, ub) domain_oracle = build_domain_oracle(Ex_mat, n) + heu = Boscia.Heuristic(Boscia.probability_rounding, 0.8, :probability_rounding) - x, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, sharpness_exponent=θ, sharpness_constant=M, domain_oracle=domain_oracle) + x, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, sharpness_exponent=θ, sharpness_constant=M, domain_oracle=domain_oracle, custom_heuristics=[heu]) + sol_a = f(x) f, grad! = build_a_criterion(Ex_mat, build_safe=false) blmo = build_blmo(m, N, ub) x0, active_set = build_start_point(Ex_mat, N, ub) z = greedy_incumbent(Ex_mat, N, ub) domain_oracle = build_domain_oracle(Ex_mat, n) + heu = Boscia.Heuristic(Boscia.probability_rounding, 0.6, :probability_rounding) - x_s, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, line_search=FrankWolfe.Secant(40, 1e-8, domain_oracle), sharpness_exponent=θ, sharpness_constant=M, domain_oracle=domain_oracle) + x_s, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, line_search=FrankWolfe.Secant(40, 1e-8, domain_oracle), sharpness_exponent=θ, sharpness_constant=M, domain_oracle=domain_oracle, custom_heuristics=[heu]) - @test isapprox(f(x), f(x_s)) + @test isapprox(sol_a, f(x_s), atol=1e-4, rtol=1e-2) end ## D-Optimal Design Problem @@ -93,9 +96,11 @@ end x0, active_set = build_start_point(Ex_mat, N, ub) z = greedy_incumbent(Ex_mat, N, ub) domain_oracle = build_domain_oracle(Ex_mat, n) + heu = Boscia.Heuristic(Boscia.probability_rounding, 0.6, :probability_rounding) - x, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, sharpness_exponent=θ, sharpness_constant=M, domain_oracle=domain_oracle) - + x, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, domain_oracle=domain_oracle, custom_heuristics=[heu]) + sol_a = f(x) #sharpness_exponent=θ, sharpness_constant=M, +@show x Ex_mat, n, N, ub = build_data(seed, m) f, grad! = build_d_criterion(Ex_mat, build_safe=false) @@ -103,8 +108,11 @@ end x0, active_set = build_start_point(Ex_mat, N, ub) z = greedy_incumbent(Ex_mat, N, ub) domain_oracle = build_domain_oracle(Ex_mat, n) + heu = Boscia.Heuristic(Boscia.probability_rounding, 0.6, :probability_rounding) - x_s, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, line_search=FrankWolfe.Secant(40, 1e-8, domain_oracle), sharpness_exponent=θ, sharpness_constant=M, domain_oracle=domain_oracle) - - @test isapprox(f(x), f(x_s)) + x_s, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, line_search=FrankWolfe.Secant(40, 1e-8, domain_oracle), domain_oracle=domain_oracle, custom_heuristics=[heu]) #sharpness_exponent=θ, sharpness_constant=M, +@show x_s +@show sum(x_s-x), findall(x-> x != 0, x-x_s) +@show sum(x_s), N, findall(x-> x > 0, x_s-ub) + @test isapprox(sol_a, f(x_s), atol=1e-4, rtol=1e-2) end From 9030728610e320f6773db33b6e06df1d3383ee49 Mon Sep 17 00:00:00 2001 From: Hendrych Date: Tue, 17 Sep 2024 11:20:57 +0200 Subject: [PATCH 34/45] More debug statements. --- examples/optimal_experiment_design.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/examples/optimal_experiment_design.jl b/examples/optimal_experiment_design.jl index 5a36e38bd..946602606 100644 --- a/examples/optimal_experiment_design.jl +++ b/examples/optimal_experiment_design.jl @@ -116,3 +116,7 @@ end @show sum(x_s), N, findall(x-> x > 0, x_s-ub) @test isapprox(sol_a, f(x_s), atol=1e-4, rtol=1e-2) end + +# in interface.jl at the polishing solution +@show x_polished + @show findall(x -> x != round.(x), x_polished[tree.root.problem.integer_variables]) From 1cc7254569ed5b221449c43ea925cd229292dc2d Mon Sep 17 00:00:00 2001 From: Hendrych Date: Tue, 17 Sep 2024 12:22:51 +0200 Subject: [PATCH 35/45] If heuristic return infeasible solution, ignore it. --- src/heuristics.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/heuristics.jl b/src/heuristics.jl index 3a5d7b3a6..7aa3f24ff 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) : true + feasible = check_feasibility ? is_linear_feasible(tree.root.problem.tlmo, x_heu) && is_integer_feasible(tree, x_heu) : false if feasible val = tree.root.problem.f(x_heu) if val < min_val @@ -166,6 +166,8 @@ function probability_rounding(tree::Bonobo.BnBTree, tlmo::Boscia.TimeTrackingLMO verbose=tree.root.options[:fwVerbose], ) + @assert sum(isapprox.(x_rounded[tlmo.blmo.int_vars], round.(x_rounded[tlmo.blmo.int_vars]))) == length(tlmo.blmo.int_vars) "$(sum(isapprox.(x_rounded[tlmo.blmo.int_vars], round.(x_rounded[tlmo.blmo.int_vars])))) == $(length(tlmo.blmo.int_vars)) $(x_rounded)" + # reset LMO to node state build_LMO(tlmo, tree.root.problem.integer_variable_bounds, original_bounds, tlmo.blmo.int_vars) From 2936851cef90817567d02fd2f24b51e60ab15360 Mon Sep 17 00:00:00 2001 From: Hendrych Date: Tue, 17 Sep 2024 13:03:22 +0200 Subject: [PATCH 36/45] Minor change. --- src/tightenings.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tightenings.jl b/src/tightenings.jl index ab6b511b6..104b3a547 100644 --- a/src/tightenings.jl +++ b/src/tightenings.jl @@ -220,7 +220,7 @@ function tightening_lowerbound(tree, node, x, lower_bound) sharpness_bound = M^(- 1 / θ) * 1/ 2 * (sqrt(bound_improvement) - M / 2 * node.dual_gap^θ)^(1 / θ) + fx - node.dual_gap @debug "Sharpness: $lower_bound -> $sharpness_bound" - @assert num_fractional == 0 || sharpness_bound > lower_bound + @assert num_fractional == 0 || sharpness_bound >= lower_bound "$(num_fractional) == 0 || $(sharpness_bound) > $(lower_bound)" end lower_bound = max(strong_convexity_bound, sharpness_bound) From 961c1a783dd7601890b619e490bccc0900ae91d8 Mon Sep 17 00:00:00 2001 From: Hendrych Date: Tue, 17 Sep 2024 13:18:26 +0200 Subject: [PATCH 37/45] Clean up example. --- examples/optimal_experiment_design.jl | 31 +++++++++++---------------- 1 file changed, 12 insertions(+), 19 deletions(-) diff --git a/examples/optimal_experiment_design.jl b/examples/optimal_experiment_design.jl index 946602606..18290f0b2 100644 --- a/examples/optimal_experiment_design.jl +++ b/examples/optimal_experiment_design.jl @@ -41,7 +41,7 @@ include("oed_utils.jl") """ seed = 1234 -m = 40 +m = 60 verbose = true ## A-Optimal Design Problem @@ -62,21 +62,20 @@ verbose = true x0, active_set = build_start_point(Ex_mat, N, ub) z = greedy_incumbent(Ex_mat, N, ub) domain_oracle = build_domain_oracle(Ex_mat, n) - heu = Boscia.Heuristic(Boscia.probability_rounding, 0.8, :probability_rounding) + heu = Boscia.Heuristic(Boscia.rounding_hyperplane_heuristic, 0.7, :hyperplane_aware_rounding) x, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, sharpness_exponent=θ, sharpness_constant=M, domain_oracle=domain_oracle, custom_heuristics=[heu]) - sol_a = f(x) f, grad! = build_a_criterion(Ex_mat, build_safe=false) blmo = build_blmo(m, N, ub) x0, active_set = build_start_point(Ex_mat, N, ub) z = greedy_incumbent(Ex_mat, N, ub) domain_oracle = build_domain_oracle(Ex_mat, n) - heu = Boscia.Heuristic(Boscia.probability_rounding, 0.6, :probability_rounding) + heu = Boscia.Heuristic(Boscia.rounding_hyperplane_heuristic, 0.7, :hyperplane_aware_rounding) - x_s, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, line_search=FrankWolfe.Secant(40, 1e-8, domain_oracle), sharpness_exponent=θ, sharpness_constant=M, domain_oracle=domain_oracle, custom_heuristics=[heu]) + x_s, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, line_search=FrankWolfe.Secant(40, 1e-8, domain_oracle), sharpness_exponent=θ, sharpness_constant=M, domain_oracle=domain_oracle, custom_heuristics=[heu]) - @test isapprox(sol_a, f(x_s), atol=1e-4, rtol=1e-2) + @test isapprox(f(x), f(x_s), atol=1e-4, rtol=1e-2) end ## D-Optimal Design Problem @@ -96,11 +95,10 @@ end x0, active_set = build_start_point(Ex_mat, N, ub) z = greedy_incumbent(Ex_mat, N, ub) domain_oracle = build_domain_oracle(Ex_mat, n) - heu = Boscia.Heuristic(Boscia.probability_rounding, 0.6, :probability_rounding) + heu = Boscia.Heuristic(Boscia.rounding_hyperplane_heuristic, 0.7, :hyperplane_aware_rounding) + + x, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, domain_oracle=domain_oracle, sharpness_exponent=θ, sharpness_constant=M, custom_heuristics=[heu]) - x, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, domain_oracle=domain_oracle, custom_heuristics=[heu]) - sol_a = f(x) #sharpness_exponent=θ, sharpness_constant=M, -@show x Ex_mat, n, N, ub = build_data(seed, m) f, grad! = build_d_criterion(Ex_mat, build_safe=false) @@ -108,15 +106,10 @@ end x0, active_set = build_start_point(Ex_mat, N, ub) z = greedy_incumbent(Ex_mat, N, ub) domain_oracle = build_domain_oracle(Ex_mat, n) - heu = Boscia.Heuristic(Boscia.probability_rounding, 0.6, :probability_rounding) + heu = Boscia.Heuristic(Boscia.rounding_hyperplane_heuristic, 0.7, :hyperplane_aware_rounding) + + x_s, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, line_search=FrankWolfe.Secant(40, 1e-8, domain_oracle), domain_oracle=domain_oracle, sharpness_exponent=θ, sharpness_constant=M, custom_heuristics=[heu]) - x_s, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, line_search=FrankWolfe.Secant(40, 1e-8, domain_oracle), domain_oracle=domain_oracle, custom_heuristics=[heu]) #sharpness_exponent=θ, sharpness_constant=M, -@show x_s -@show sum(x_s-x), findall(x-> x != 0, x-x_s) -@show sum(x_s), N, findall(x-> x > 0, x_s-ub) - @test isapprox(sol_a, f(x_s), atol=1e-4, rtol=1e-2) + @test isapprox(f(x), f(x_s), atol=1e-4, rtol=1e-2) end -# in interface.jl at the polishing solution -@show x_polished - @show findall(x -> x != round.(x), x_polished[tree.root.problem.integer_variables]) From e1d298b8741779d52315ae52f91dba478e93ea90 Mon Sep 17 00:00:00 2001 From: Hendrych Date: Tue, 17 Sep 2024 13:41:18 +0200 Subject: [PATCH 38/45] Show solutions. --- examples/optimal_experiment_design.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/examples/optimal_experiment_design.jl b/examples/optimal_experiment_design.jl index 18290f0b2..027b4f45e 100644 --- a/examples/optimal_experiment_design.jl +++ b/examples/optimal_experiment_design.jl @@ -75,6 +75,8 @@ verbose = true x_s, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, line_search=FrankWolfe.Secant(40, 1e-8, domain_oracle), sharpness_exponent=θ, sharpness_constant=M, domain_oracle=domain_oracle, custom_heuristics=[heu]) + @show x + @show x_s @test isapprox(f(x), f(x_s), atol=1e-4, rtol=1e-2) end @@ -99,7 +101,7 @@ end x, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, domain_oracle=domain_oracle, sharpness_exponent=θ, sharpness_constant=M, custom_heuristics=[heu]) - Ex_mat, n, N, ub = build_data(seed, m) + # Ex_mat, n, N, ub = build_data(seed, m) f, grad! = build_d_criterion(Ex_mat, build_safe=false) blmo = build_blmo(m, N, ub) @@ -110,6 +112,8 @@ end x_s, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, line_search=FrankWolfe.Secant(40, 1e-8, domain_oracle), domain_oracle=domain_oracle, sharpness_exponent=θ, sharpness_constant=M, custom_heuristics=[heu]) + @show x + @show x_s @test isapprox(f(x), f(x_s), atol=1e-4, rtol=1e-2) end From 4e42b57f7d2fd48eef9e46c64a6c14d476a373bd Mon Sep 17 00:00:00 2001 From: Hendrych Date: Tue, 17 Sep 2024 14:21:58 +0200 Subject: [PATCH 39/45] Don't use vertex storage for the heuristics. --- src/heuristics.jl | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/src/heuristics.jl b/src/heuristics.jl index 7aa3f24ff..c30392a8c 100644 --- a/src/heuristics.jl +++ b/src/heuristics.jl @@ -159,14 +159,11 @@ function probability_rounding(tree::Bonobo.BnBTree, tlmo::Boscia.TimeTrackingLMO line_search=tree.root.options[:lineSearch], lazy=tree.root.options[:lazy], lazy_tolerance=tree.root.options[:lazy_tolerance], - add_dropped_vertices=tree.root.options[:use_shadow_set], - use_extra_vertex_storage=tree.root.options[:use_shadow_set], - extra_vertex_storage=node.discarded_vertices, callback=tree.root.options[:callback], verbose=tree.root.options[:fwVerbose], ) - @assert sum(isapprox.(x_rounded[tlmo.blmo.int_vars], round.(x_rounded[tlmo.blmo.int_vars]))) == length(tlmo.blmo.int_vars) "$(sum(isapprox.(x_rounded[tlmo.blmo.int_vars], round.(x_rounded[tlmo.blmo.int_vars])))) == $(length(tlmo.blmo.int_vars)) $(x_rounded)" + @assert sum(isapprox.(x_rounded[tlmo.blmo.int_vars], round.(x_rounded[tlmo.blmo.int_vars]))) == length(tlmo.blmo.int_vars) "$(sum(isapprox.(x_rounded[tlmo.blmo.int_vars], round.(x_rounded[tlmo.blmo.int_vars])))) == $(length(tlmo.blmo.int_vars)) $(x_rounded[tlmo.blmo.int_vars])" # reset LMO to node state build_LMO(tlmo, tree.root.problem.integer_variable_bounds, original_bounds, tlmo.blmo.int_vars) From 98378586b1621ab0fc95aed43af209f5fd292e6c Mon Sep 17 00:00:00 2001 From: Hendrych Date: Wed, 18 Sep 2024 16:06:59 +0200 Subject: [PATCH 40/45] Final clean up. --- examples/oed_utils.jl | 4 +-- examples/optimal_experiment_design.jl | 36 +++++++++++++++------------ 2 files changed, 21 insertions(+), 19 deletions(-) diff --git a/examples/oed_utils.jl b/examples/oed_utils.jl index b67bfb42b..32064d071 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(seed, m) +function build_data(m) # set up - Random.seed!(seed) - n = Int(floor(m/10)) B = rand(m,n) diff --git a/examples/optimal_experiment_design.jl b/examples/optimal_experiment_design.jl index 027b4f45e..ca6ebfbba 100644 --- a/examples/optimal_experiment_design.jl +++ b/examples/optimal_experiment_design.jl @@ -40,7 +40,6 @@ include("oed_utils.jl") https://github.com/ZIB-IOL/FrankWolfe.jl/blob/master/examples/optimal_experiment_design.jl """ -seed = 1234 m = 60 verbose = true @@ -48,7 +47,7 @@ verbose = true @testset "A-Optimal Design" begin - Ex_mat, n, N, ub = build_data(seed, m) + Ex_mat, n, N, ub = build_data(m) # sharpness constants σ = minimum(Ex_mat' * Ex_mat) @@ -57,33 +56,40 @@ verbose = true M = n * σ^4 / λ_max^2 - f, grad! = build_a_criterion(Ex_mat, build_safe=true) + g, grad! = build_a_criterion(Ex_mat, build_safe=true) blmo = build_blmo(m, N, ub) x0, active_set = build_start_point(Ex_mat, N, ub) z = greedy_incumbent(Ex_mat, N, ub) domain_oracle = build_domain_oracle(Ex_mat, n) heu = Boscia.Heuristic(Boscia.rounding_hyperplane_heuristic, 0.7, :hyperplane_aware_rounding) - x, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, sharpness_exponent=θ, sharpness_constant=M, domain_oracle=domain_oracle, custom_heuristics=[heu]) + x, _, result = Boscia.solve(g, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, domain_oracle=domain_oracle, sharpness_exponent=θ, sharpness_constant=M, custom_heuristics=[heu]) #sharpness_exponent=θ, sharpness_constant=M, + + gradient = similar(x) + grad!(gradient, x) + v = Boscia.compute_extreme_point(blmo, gradient) + dual_gap = FrankWolfe.dot(gradient, x) - FrankWolfe.dot(gradient, v) + @show dual_gap - f, grad! = build_a_criterion(Ex_mat, build_safe=false) blmo = build_blmo(m, N, ub) x0, active_set = build_start_point(Ex_mat, N, ub) z = greedy_incumbent(Ex_mat, N, ub) domain_oracle = build_domain_oracle(Ex_mat, n) heu = Boscia.Heuristic(Boscia.rounding_hyperplane_heuristic, 0.7, :hyperplane_aware_rounding) - x_s, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, line_search=FrankWolfe.Secant(40, 1e-8, domain_oracle), sharpness_exponent=θ, sharpness_constant=M, domain_oracle=domain_oracle, custom_heuristics=[heu]) + x_s, _, result = Boscia.solve(g, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, line_search=FrankWolfe.Secant(40, 1e-8, domain_oracle), sharpness_exponent=θ, sharpness_constant=M, domain_oracle=domain_oracle, custom_heuristics=[heu]) #sharpness_exponent=θ, sharpness_constant=M, @show x @show x_s - @test isapprox(f(x), f(x_s), atol=1e-4, rtol=1e-2) -end + #@show f(x), f(x_s) + @show g(x), g(x_s) + @test isapprox(g(x), g(x_s), atol=1e-4, rtol=5e-2) +end ## D-Optimal Design Problem @testset "D-optimal Design" begin - Ex_mat, n, N, ub = build_data(seed, m) + Ex_mat, n, N, ub = build_data(m) # sharpness constants σ = minimum(Ex_mat' * Ex_mat) @@ -92,28 +98,26 @@ end M = n * σ^4 / (2 * λ_max^2) - f, grad! = build_d_criterion(Ex_mat, build_safe=true) + g, grad! = build_d_criterion(Ex_mat, build_safe=true) blmo = build_blmo(m, N, ub) x0, active_set = build_start_point(Ex_mat, N, ub) z = greedy_incumbent(Ex_mat, N, ub) domain_oracle = build_domain_oracle(Ex_mat, n) heu = Boscia.Heuristic(Boscia.rounding_hyperplane_heuristic, 0.7, :hyperplane_aware_rounding) - x, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, domain_oracle=domain_oracle, sharpness_exponent=θ, sharpness_constant=M, custom_heuristics=[heu]) - - # Ex_mat, n, N, ub = build_data(seed, m) + x, _, result = Boscia.solve(g, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, domain_oracle=domain_oracle, sharpness_exponent=θ, sharpness_constant=M, custom_heuristics=[heu]) - f, grad! = build_d_criterion(Ex_mat, build_safe=false) blmo = build_blmo(m, N, ub) x0, active_set = build_start_point(Ex_mat, N, ub) z = greedy_incumbent(Ex_mat, N, ub) domain_oracle = build_domain_oracle(Ex_mat, n) heu = Boscia.Heuristic(Boscia.rounding_hyperplane_heuristic, 0.7, :hyperplane_aware_rounding) - x_s, _, result = Boscia.solve(f, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, line_search=FrankWolfe.Secant(40, 1e-8, domain_oracle), domain_oracle=domain_oracle, sharpness_exponent=θ, sharpness_constant=M, custom_heuristics=[heu]) + x_s, _, result = Boscia.solve(g, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, line_search=FrankWolfe.Secant(40, 1e-8, domain_oracle), domain_oracle=domain_oracle, sharpness_exponent=θ, sharpness_constant=M, custom_heuristics=[heu]) @show x @show x_s - @test isapprox(f(x), f(x_s), atol=1e-4, rtol=1e-2) + @show g(x), g(x_s) + @test isapprox(g(x), g(x_s), atol=1e-4, rtol=1e-2) end From c0c822fa5f6a80151f1dbf4e8ddfc31fef494c37 Mon Sep 17 00:00:00 2001 From: Hendrych Date: Wed, 18 Sep 2024 17:39:55 +0200 Subject: [PATCH 41/45] Run tests with one specific seed. --- examples/birkhoff.jl | 6 ------ examples/int_sparse_reg.jl | 6 ------ examples/mps-example.jl | 4 ---- examples/mps-examples/mip-examples.jl | 4 ---- examples/nonlinear.jl | 4 ---- examples/poisson_reg.jl | 6 ------ examples/portfolio.jl | 4 ---- 7 files changed, 34 deletions(-) diff --git a/examples/birkhoff.jl b/examples/birkhoff.jl index 7c7025b48..655ef7ccf 100644 --- a/examples/birkhoff.jl +++ b/examples/birkhoff.jl @@ -12,12 +12,6 @@ import HiGHS # https://arxiv.org/pdf/2011.02752.pdf # https://www.sciencedirect.com/science/article/pii/S0024379516001257 -# For bug hunting: -seed = rand(UInt64) -@show seed -seed = 0x3eb09305cecf69f0 -Random.seed!(seed) - # min_{X, θ} 1/2 * || ∑_{i in [k]} θ_i X_i - Xhat ||^2 # θ ∈ Δ_k (simplex) diff --git a/examples/int_sparse_reg.jl b/examples/int_sparse_reg.jl index 1d7a0a3b9..9d6a037d8 100644 --- a/examples/int_sparse_reg.jl +++ b/examples/int_sparse_reg.jl @@ -21,12 +21,6 @@ const MOI = MathOptInterface # r - controls how often we have to maximal split on a index. # k - is the sparsity parameter. We only want a few non zero entries. -# For bug hunting: -seed = rand(UInt64) -@show seed -#seed = 0xeadb922ca734998b -Random.seed!(seed) - n = 10 m = 30 l = 5 diff --git a/examples/mps-example.jl b/examples/mps-example.jl index 6ea84f0ed..e60208919 100644 --- a/examples/mps-example.jl +++ b/examples/mps-example.jl @@ -16,10 +16,6 @@ const MOI = MathOptInterface # Number of binaries 231 # Number of constraints 198 -seed = rand(UInt64) -@show seed -Random.seed!(seed) - src = MOI.FileFormats.Model(filename="22433.mps") MOI.read_from_file(src, joinpath(@__DIR__, "mps-examples/mps-files/22433.mps")) diff --git a/examples/mps-examples/mip-examples.jl b/examples/mps-examples/mip-examples.jl index 173dfd089..aba7854be 100644 --- a/examples/mps-examples/mip-examples.jl +++ b/examples/mps-examples/mip-examples.jl @@ -21,10 +21,6 @@ import Ipopt # ran14x18-disj-8 https://miplib.zib.de/instance_details_ran14x18-disj-8.html # timtab1 https://miplib.zib.de/instance_details_timtab1.html (Takes LONG!) -seed = rand(UInt64) -@show seed -Random.seed!(seed) - # To see debug statements #ENV["JULIA_DEBUG"] = "Boscia" diff --git a/examples/nonlinear.jl b/examples/nonlinear.jl index 4fac4e8f9..9ceb4b45f 100644 --- a/examples/nonlinear.jl +++ b/examples/nonlinear.jl @@ -11,10 +11,6 @@ using Dates # using SCIP # const MOI = MathOptInterface -n = 40 -seed = 10 - -Random.seed!(seed) ################################################################ # alternative implementation of LMO using MOI and SCIP diff --git a/examples/poisson_reg.jl b/examples/poisson_reg.jl index 1f87527f1..8bde3f90e 100644 --- a/examples/poisson_reg.jl +++ b/examples/poisson_reg.jl @@ -11,12 +11,6 @@ const MOI = MathOptInterface # Poisson sparse regression -# For bug hunting: -seed = rand(UInt64) -@show seed -#seed = 0xfe03ee83ca373eab -Random.seed!(seed) - # min_{w, b, z} ∑_i exp(w x_i + b) - y_i (w x_i + b) + α norm(w)^2 # s.t. -N z_i <= w_i <= N z_i # b ∈ [-N, N] diff --git a/examples/portfolio.jl b/examples/portfolio.jl index 8b0221384..22b4c3c22 100644 --- a/examples/portfolio.jl +++ b/examples/portfolio.jl @@ -7,10 +7,6 @@ using LinearAlgebra import MathOptInterface const MOI = MathOptInterface -# seed = 0x946d4b7835e92ffa takes 90 minutes to solve! -> not anymore -seed = 0x946d4b7835e92ffa -Random.seed!(seed) - n = 30 const ri = rand(n) const ai = rand(n) From 7d93f84b2b30aceb2ba69e352a8112227ae2b002 Mon Sep 17 00:00:00 2001 From: Hendrych Date: Wed, 18 Sep 2024 18:33:21 +0200 Subject: [PATCH 42/45] Added time limit. --- examples/nonlinear.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/nonlinear.jl b/examples/nonlinear.jl index 9ceb4b45f..fdc27539d 100644 --- a/examples/nonlinear.jl +++ b/examples/nonlinear.jl @@ -109,6 +109,6 @@ heu2 = Boscia.Heuristic(Boscia.rounding_lmo_01_heuristic, 0.8, :lmo_rounding) heuristics = [heu, heu2] # heuristics = [] -x, _, _ = Boscia.solve(f, grad!, lmo, verbose=true, print_iter=500, custom_heuristics=heuristics) +x, _, _ = Boscia.solve(f, grad!, lmo, verbose=true, print_iter=500, custom_heuristics=heuristics, time_limit=300) @show x From c5b00196600c3a7b3228b6c28b7b111dcd8f1b7c Mon Sep 17 00:00:00 2001 From: Hendrych Date: Thu, 19 Sep 2024 17:30:44 +0200 Subject: [PATCH 43/45] Corrected sharpness constant. --- examples/optimal_experiment_design.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/optimal_experiment_design.jl b/examples/optimal_experiment_design.jl index ca6ebfbba..a959188d5 100644 --- a/examples/optimal_experiment_design.jl +++ b/examples/optimal_experiment_design.jl @@ -53,7 +53,7 @@ verbose = true σ = minimum(Ex_mat' * Ex_mat) λ_max = maximum(ub) * maximum([norm(Ex_mat[i,:])^2 for i=1:size(Ex_mat,1)]) θ = 1/2 - M = n * σ^4 / λ_max^2 + M = n * σ^4 / λ_max^3 g, grad! = build_a_criterion(Ex_mat, build_safe=true) @@ -95,7 +95,7 @@ end σ = minimum(Ex_mat' * Ex_mat) λ_max = maximum(ub) * maximum([norm(Ex_mat[i,:])^2 for i=1:size(Ex_mat,1)]) θ = 1/2 - M = n * σ^4 / (2 * λ_max^2) + M = n * σ^4 / λ_max^2 g, grad! = build_d_criterion(Ex_mat, build_safe=true) From 069c4d83d5ef0310267f7545d098fe2da9d8278b Mon Sep 17 00:00:00 2001 From: Hendrych Date: Thu, 19 Sep 2024 18:04:30 +0200 Subject: [PATCH 44/45] Corrected sharpness constant and slight relaxation of test for A-Optimal. --- examples/optimal_experiment_design.jl | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/examples/optimal_experiment_design.jl b/examples/optimal_experiment_design.jl index a959188d5..236c40939 100644 --- a/examples/optimal_experiment_design.jl +++ b/examples/optimal_experiment_design.jl @@ -39,8 +39,11 @@ include("oed_utils.jl") A continuous version of the problem can be found in the examples in FrankWolfe.jl: https://github.com/ZIB-IOL/FrankWolfe.jl/blob/master/examples/optimal_experiment_design.jl """ +seed = rand(UInt64) -m = 60 +@show seed +Random.seed!(seed) +m = 40 verbose = true ## A-Optimal Design Problem @@ -49,11 +52,13 @@ verbose = true Ex_mat, n, N, ub = build_data(m) + @show Ex_mat + # sharpness constants σ = minimum(Ex_mat' * Ex_mat) λ_max = maximum(ub) * maximum([norm(Ex_mat[i,:])^2 for i=1:size(Ex_mat,1)]) θ = 1/2 - M = n * σ^4 / λ_max^3 + M = sqrt(λ_max^3/ n * σ^4) g, grad! = build_a_criterion(Ex_mat, build_safe=true) @@ -61,9 +66,10 @@ verbose = true x0, active_set = build_start_point(Ex_mat, N, ub) z = greedy_incumbent(Ex_mat, N, ub) domain_oracle = build_domain_oracle(Ex_mat, n) + line_search = FrankWolfe.MonotonicGenericStepsize(FrankWolfe.Adaptive(), domain_oracle) heu = Boscia.Heuristic(Boscia.rounding_hyperplane_heuristic, 0.7, :hyperplane_aware_rounding) - x, _, result = Boscia.solve(g, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, domain_oracle=domain_oracle, sharpness_exponent=θ, sharpness_constant=M, custom_heuristics=[heu]) #sharpness_exponent=θ, sharpness_constant=M, + x, _, result = Boscia.solve(g, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, domain_oracle=domain_oracle, custom_heuristics=[heu], sharpness_exponent=θ, sharpness_constant=M, line_search=line_search) #sharpness_exponent=θ, sharpness_constant=M, gradient = similar(x) grad!(gradient, x) @@ -77,13 +83,13 @@ verbose = true domain_oracle = build_domain_oracle(Ex_mat, n) heu = Boscia.Heuristic(Boscia.rounding_hyperplane_heuristic, 0.7, :hyperplane_aware_rounding) - x_s, _, result = Boscia.solve(g, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, line_search=FrankWolfe.Secant(40, 1e-8, domain_oracle), sharpness_exponent=θ, sharpness_constant=M, domain_oracle=domain_oracle, custom_heuristics=[heu]) #sharpness_exponent=θ, sharpness_constant=M, + x_s, _, result = Boscia.solve(g, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, line_search=FrankWolfe.Secant(40, 1e-8, domain_oracle), domain_oracle=domain_oracle, sharpness_exponent=θ, sharpness_constant=M, custom_heuristics=[heu]) #sharpness_exponent=θ, sharpness_constant=M, @show x @show x_s #@show f(x), f(x_s) @show g(x), g(x_s) - @test isapprox(g(x), g(x_s), atol=1e-4, rtol=5e-2) + @test isapprox(g(x), g(x_s), atol=1e-3, rtol=5e-2) end ## D-Optimal Design Problem @@ -95,7 +101,7 @@ end σ = minimum(Ex_mat' * Ex_mat) λ_max = maximum(ub) * maximum([norm(Ex_mat[i,:])^2 for i=1:size(Ex_mat,1)]) θ = 1/2 - M = n * σ^4 / λ_max^2 + M = sqrt(2 * λ_max^2/ n * σ^4 ) g, grad! = build_d_criterion(Ex_mat, build_safe=true) From db20f7e68fc907fd1c2cbe51082a309c2e1088db Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Fri, 20 Sep 2024 12:01:35 +0200 Subject: [PATCH 45/45] Update src/tightenings.jl --- src/tightenings.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/tightenings.jl b/src/tightenings.jl index 104b3a547..2cde12a4f 100644 --- a/src/tightenings.jl +++ b/src/tightenings.jl @@ -229,7 +229,8 @@ function tightening_lowerbound(tree, node, x, lower_bound) end """ -Use strong convexity and/or sharpness to potentially remove one of the children nodes +Use strong convexity and/or sharpness to potentially remove one of the children nodes. +If both sharpness and strong convexity parameters are provided, strong convexity is preferred. """ function prune_children(tree, node, lower_bound_base, x, vidx) prune_left = false