diff --git a/examples/approx_planted_point.jl b/examples/approx_planted_point.jl index d027eab42..04a709ea0 100644 --- a/examples/approx_planted_point.jl +++ b/examples/approx_planted_point.jl @@ -9,6 +9,8 @@ using Distributions import MathOptInterface const MOI = MathOptInterface +include("cube_blmo.jl") + n = 20 diffi = Random.rand(Bool, n) * 0.6 .+ 0.3 @@ -47,7 +49,7 @@ diffi = Random.rand(Bool, n) * 0.6 .+ 0.3 push!(bounds, (i, 0.0), :greaterthan) push!(bounds, (i, 1.0), :lessthan) end - blmo = Boscia.CubeBLMO(n, int_vars, bounds) + blmo = CubeBLMO(n, int_vars, bounds) x, _, result = Boscia.solve(f, grad!, blmo, verbose=true) @@ -110,7 +112,7 @@ end push!(bounds, (i, 0.0), :greaterthan) push!(bounds, (i, 1.0), :lessthan) end - blmo = Boscia.CubeBLMO(n, int_vars, bounds) + blmo = CubeBLMO(n, int_vars, bounds) x, _, result = Boscia.solve(f, grad!, blmo, verbose=true) diff --git a/src/cube_blmo.jl b/examples/cube_blmo.jl similarity index 54% rename from src/cube_blmo.jl rename to examples/cube_blmo.jl index 9ac7f8fec..76fa6230d 100644 --- a/src/cube_blmo.jl +++ b/examples/cube_blmo.jl @@ -1,13 +1,17 @@ +## How to implement the BLMO Interface using the cube as an example +using Boscia +using Bonobo +using Dates """ CubeBLMO A Bounded Linear Minimization Oracle over a cube. """ -mutable struct CubeBLMO <: BoundedLinearMinimizationOracle +mutable struct CubeBLMO <: Boscia.BoundedLinearMinimizationOracle n::Int int_vars::Vector{Int} - bounds::IntegerBounds + bounds::Boscia.IntegerBounds solving_time::Float64 end @@ -16,7 +20,7 @@ CubeBLMO(n, int_vars, bounds) = CubeBLMO(n, int_vars, bounds, 0.0) ## Necessary # computing an extreme point for the cube amounts to checking the sign of the gradient -function compute_extreme_point(blmo::CubeBLMO, d; kwargs...) +function Boscia.compute_extreme_point(blmo::CubeBLMO, d; kwargs...) time_ref = Dates.now() v = zeros(length(d)) for i in eachindex(d) @@ -28,8 +32,8 @@ end ## -function build_global_bounds(blmo::CubeBLMO, integer_variables) - global_bounds = IntegerBounds() +function Boscia.build_global_bounds(blmo::CubeBLMO, integer_variables) + global_bounds = Boscia.IntegerBounds() for i in 1:blmo.n if i in integer_variables push!(global_bounds, (i, blmo.bounds[i, :lessthan]), :lessthan) @@ -41,42 +45,34 @@ end ## Read information from problem -function get_list_of_variables(blmo::CubeBLMO) +function Boscia.get_list_of_variables(blmo::CubeBLMO) return blmo.n, collect(1:blmo.n) end # Get list of integer variables, respectively. -function get_integer_variables(blmo::CubeBLMO) +function Boscia.get_integer_variables(blmo::CubeBLMO) return blmo.int_vars end -function get_int_var(blmo::CubeBLMO, cidx) +function Boscia.get_int_var(blmo::CubeBLMO, cidx) return cidx end -function get_lower_bound_list(blmo::CubeBLMO) +function Boscia.get_lower_bound_list(blmo::CubeBLMO) return keys(blmo.bounds.lower_bounds) end -function get_upper_bound_list(blmo::CubeBLMO) +function Boscia.get_upper_bound_list(blmo::CubeBLMO) return keys(blmo.bounds.upper_bounds) end -function get_bound(blmo::CubeBLMO, c_idx, sense::Symbol) +function Boscia.get_bound(blmo::CubeBLMO, c_idx, sense::Symbol) @assert sense == :lessthan || sense == :greaterthan return blmo[c_idx, sense] end -#function get_lower_bound(blmo::CubeBLMO, c_idx) -# return blmo.bounds[c_idx, :greaterthan] -#end - -#function get_upper_bound(blmo::CubeBLMO, c_idx) -# return blmo.bounds[c_idx, :lessthan] -#end - ## Changing the bounds constraints. -function set_bound!(blmo::CubeBLMO, c_idx, value, sense::Symbol) +function Boscia.set_bound!(blmo::CubeBLMO, c_idx, value, sense::Symbol) if sense == :greaterthan blmo.bounds.lower_bounds[c_idx] = value elseif sense == :lessthan @@ -86,29 +82,29 @@ function set_bound!(blmo::CubeBLMO, c_idx, value, sense::Symbol) end end -function delete_bounds!(blmo::CubeBLMO, cons_delete) +function Boscia.delete_bounds!(blmo::CubeBLMO, cons_delete) # For the cube this shouldn't happen! Otherwise we get unbounded! if !isempty(cons_delete) error("Trying to delete bounds of the cube!") end end -function add_bound_constraint!(blmo::CubeBLMO, key, value, sense::Symbol) +function Boscia.add_bound_constraint!(blmo::CubeBLMO, key, value, sense::Symbol) # Should not be necessary return error("Trying to add bound constraints of the cube!") end ## Checks -function is_constraint_on_int_var(blmo::CubeBLMO, c_idx, int_vars) +function Boscia.is_constraint_on_int_var(blmo::CubeBLMO, c_idx, int_vars) return c_idx in int_vars end -function is_bound_in(blmo::CubeBLMO, c_idx, bounds) +function Boscia.is_bound_in(blmo::CubeBLMO, c_idx, bounds) return haskey(bounds, c_idx) end -function is_linear_feasible(blmo::CubeBLMO, v::AbstractVector) +function Boscia.is_linear_feasible(blmo::CubeBLMO, v::AbstractVector) for i in eachindex(v) if !( blmo.bounds[i, :greaterthan] ≤ v[i] + 1e-6 || @@ -123,16 +119,15 @@ function is_linear_feasible(blmo::CubeBLMO, v::AbstractVector) return true end -function has_integer_constraint(blmo::CubeBLMO, idx) +function Boscia.has_integer_constraint(blmo::CubeBLMO, idx) return idx in blmo.int_vars end - ###################### Optional ## Safety Functions -function build_LMO_correct(blmo::CubeBLMO, node_bounds) +function Boscia.build_LMO_correct(blmo::CubeBLMO, node_bounds) for key in keys(node_bounds.lower_bounds) if !haskey(blmo.bounds, (key, :greaterthan)) || blmo.bounds[key, :greaterthan] != node_bounds[key, :greaterthan] @@ -148,60 +143,23 @@ function build_LMO_correct(blmo::CubeBLMO, node_bounds) return true end -function check_feasibility(blmo::CubeBLMO) +function Boscia.check_feasibility(blmo::CubeBLMO) for i in 1:blmo.n if !haskey(blmo.bounds, (i, :greaterthan)) || !haskey(blmo.bounds, (i, :lessthan)) - return UNBOUNDED + return Boscia.UNBOUNDED end if blmo.bounds[i, :greaterthan] > blmo.bounds[i, :lessthan] - return INFEASIBLE + return Boscia.INFEASIBLE end end - return OPTIMAL + return Boscia.OPTIMAL end -function is_valid_split(tree::Bonobo.BnBTree, blmo::CubeBLMO, vidx::Int) +function Boscia.is_valid_split(tree::Bonobo.BnBTree, blmo::CubeBLMO, vidx::Int) return blmo.bounds[vidx, :lessthan] != blmo.bounds[vidx, :greaterthan] end ## Logs -function get_BLMO_solve_data(blmo::CubeBLMO) +function Boscia.get_BLMO_solve_data(blmo::CubeBLMO) return blmo.solving_time, 0.0, 0.0 -end - -######################################################################## -""" - CubeSimpleBLMO{T}(lower_bounds, upper_bounds) - -Hypercube with lower and upper bounds implementing the `SimpleBoundableLMO` interface. -""" -struct CubeSimpleBLMO <: SimpleBoundableLMO - lower_bounds::Vector{Float64} - upper_bounds::Vector{Float64} - int_vars::Vector{Int} -end - -function bounded_compute_extreme_point(sblmo::CubeSimpleBLMO, d, lb, ub, int_vars; kwargs...) - v = zeros(length(d)) - for i in eachindex(d) - if i in int_vars - idx = findfirst(x -> x == i, int_vars) - v[i] = d[i] > 0 ? lb[idx] : ub[idx] - else - v[i] = d[i] > 0 ? sblmo.lower_bounds[i] : sblmo.upper_bounds[i] - end - end - return v -end - -function is_linear_feasible(sblmo::CubeSimpleBLMO, v) - for i in setdiff(eachindex(v), sblmo.int_vars) - if !(sblmo.lower_bounds[i] ≤ v[i] + 1e-6 || !(v[i] - 1e-6 ≤ blmo.upper_bounds[i])) - @debug( - "Vertex entry: $(v[i]) Lower bound: $(blmo.bounds[i, :greaterthan]) Upper bound: $(blmo.bounds[i, :lessthan]))" - ) - return false - end - end - return true -end +end \ No newline at end of file diff --git a/src/Boscia.jl b/src/Boscia.jl index 3b5a14cf0..5e92fc6ca 100644 --- a/src/Boscia.jl +++ b/src/Boscia.jl @@ -27,7 +27,7 @@ include("utilities.jl") include("interface.jl") include("managed_blmo.jl") include("MOI_bounded_oracle.jl") -include("cube_blmo.jl") +include("polytope_blmos.jl") # For extensions if !isdefined(Base, :get_extension) diff --git a/src/managed_blmo.jl b/src/managed_blmo.jl index 6c279e490..37c9f4f58 100644 --- a/src/managed_blmo.jl +++ b/src/managed_blmo.jl @@ -174,7 +174,7 @@ function is_linear_feasible(blmo::ManagedBoundedLMO, v::AbstractVector) blmo.lower_bounds[i] ≤ v[int_var] + 1e-6 || !(v[int_var] - 1e-6 ≤ blmo.upper_bounds[i]) ) @debug( - "Vertex entry: $(v[int_var]) Lower bound: $(blmo.lower_bounds[i]) Upper bound: $(blmo.upper_bounds[i]))" + "Variable: $(int_var) Vertex entry: $(v[int_var]) Lower bound: $(blmo.lower_bounds[i]) Upper bound: $(blmo.upper_bounds[i]))" ) return false end diff --git a/src/polytope_blmos.jl b/src/polytope_blmos.jl new file mode 100644 index 000000000..4bda95e36 --- /dev/null +++ b/src/polytope_blmos.jl @@ -0,0 +1,120 @@ +""" + CubeSimpleBLMO{T}(lower_bounds, upper_bounds) + +Hypercube with lower and upper bounds implementing the `SimpleBoundableLMO` interface. +""" +struct CubeSimpleBLMO <: SimpleBoundableLMO + lower_bounds::Vector{Float64} + upper_bounds::Vector{Float64} + int_vars::Vector{Int} +end + +""" +If the entry is positve, choose the lower bound. Else, choose the upper bound. +""" +function bounded_compute_extreme_point(sblmo::CubeSimpleBLMO, d, lb, ub, int_vars; kwargs...) + v = zeros(length(d)) + for i in eachindex(d) + if i in int_vars + idx = findfirst(x -> x == i, int_vars) + v[i] = d[i] > 0 ? lb[idx] : ub[idx] + else + v[i] = d[i] > 0 ? sblmo.lower_bounds[i] : sblmo.upper_bounds[i] + end + end + return v +end + +function is_linear_feasible(sblmo::CubeSimpleBLMO, v) + for i in setdiff(eachindex(v), sblmo.int_vars) + if !(sblmo.lower_bounds[i] ≤ v[i] + 1e-6 || !(v[i] - 1e-6 ≤ blmo.upper_bounds[i])) + @debug( + "Vertex entry: $(v[i]) Lower bound: $(blmo.bounds[i, :greaterthan]) Upper bound: $(blmo.bounds[i, :lessthan]))" + ) + return false + end + end + return true +end + +""" + ProbablitySimplexSimpleBLMO(N) + +Scaled Probability Simplex: ∑ x = N. +""" +struct ProbabilitySimplexSimpleBLMO <: SimpleBoundableLMO + N::Float64 +end + +""" +Assign the largest possible values to the entries corresponding to the smallest entries of d. +""" +function bounded_compute_extreme_point(sblmo::ProbabilitySimplexSimpleBLMO, d, lb, ub, int_vars; kwargs...) + v = zeros(length(d)) + indices = collect(1:length(d)) + perm = sortperm(d) + + # The lower bounds always have to be met. + v[int_vars] = lb + + for i in indices[perm] + if i in int_vars + idx = findfirst(x -> x == i, int_vars) + v[i] += min(ub[idx]-lb[idx], sblmo.N - sum(v)) + else + v[i] += N - sum(v) + end + end + return v +end + +function is_linear_feasible(sblmo::ProbabilitySimplexSimpleBLMO, v) + if sum(v .≥ 0) < length(v) + @debug "v has negative entries: $(v)" + return false + end + return isapprox(sum(v), sblmo.N, atol=1e-4, rtol=1e-2) +end + +""" + UnitSimplexSimpleBLMO(N) + +Scaled Unit Simplex: ∑ x ≤ N. +""" +struct UnitSimplexSimpleBLMO <: SimpleBoundableLMO + N::Float64 +end + +""" +For all positive entries of d, assign the corresponding lower bound. +For non-positive entries, assign largest possible value in increasing order. +""" +function bounded_compute_extreme_point(sblmo::UnitSimplexSimpleBLMO, d, lb, ub, int_vars; kwargs...) + v = zeros(length(d)) + # The wloer bounds always have to be met. + v[int_vars] = lb + cont_vars = setdiff(collect(1:length(d)), int_vars) + if !isempty(cont_vars) + v[cont_vars] .= 0.0 + end + + idx_neg = findall(x-> x <= 0, d) + perm = sortperm(d[idx_neg]) + for i in idx_neg[perm] + if i in int_vars + idx = findfirst(x -> x == i, int_vars) + v[i] += min(ub[idx]-lb[idx], sblmo.N - sum(v)) + else + v[i] += N - sum(v) + end + end + return v +end + +function is_linear_feasible(sblmo::UnitSimplexSimpleBLMO, v) + if sum(v .≥ 0) < length(v) + @debug "v has negative entries: $(v)" + return false + end + return sum(v) ≤ sblmo.N + 1e-3 +end diff --git a/test/LMO_test.jl b/test/LMO_test.jl index 35005e0fd..c879bd715 100644 --- a/test/LMO_test.jl +++ b/test/LMO_test.jl @@ -11,8 +11,6 @@ using Dates const MOI = MathOptInterface const MOIU = MOI.Utilities -import MathOptSetDistances -const MOD = MathOptSetDistances @testset "Integer bounds" begin n = 10 integer_bound = Boscia.IntegerBounds() @@ -66,6 +64,27 @@ diffi = Random.rand(Bool, n) * 0.6 .+ 0.3 function grad!(storage, x) @. storage = x - diffi end + + lbs = zeros(n) + ubs = ones(n) + int_vars = collect(1:n) + + sblmo = Boscia.CubeSimpleBLMO(lbs, ubs, int_vars) + + x, _, result = + Boscia.solve(f, grad!, sblmo, lbs[int_vars], ubs[int_vars], int_vars, n) + + @test sum(isapprox.(x, round.(diffi), atol=1e-6, rtol=1e-2)) == n + @test isapprox(f(x), f(result[:raw_solution]), atol=1e-6, rtol=1e-3) +end + +@testset "BLMO - Strong Branching" 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 "Partial Strong Branching" begin int_vars = collect(1:n) lbs = zeros(n) @@ -102,3 +121,48 @@ diffi = Random.rand(Bool, n) * 0.6 .+ 0.3 @test isapprox(f(x), f(result[:raw_solution]), atol=1e-6, rtol=1e-3) end end + +n = 20 +x_sol = rand(1:Int(floor(n/4)), n) +N = sum(x_sol) +dir = vcat(fill(1, Int(floor(n/2))), fill(-1, Int(floor(n/2))), fill(0, mod(n,2))) +diffi = x_sol + 0.3 * dir + +@testset "Probability Simplex LMO" 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 + + sblmo = Boscia.ProbabilitySimplexSimpleBLMO(N) + + x, _, result = + Boscia.solve(f, grad!, sblmo, fill(0.0, n), fill(1.0*N, n), collect(1:n), n) + + @test sum(isapprox.(x, x_sol, atol=1e-6, rtol=1e-2)) == n + @test isapprox(f(x), f(result[:raw_solution]), atol=1e-6, rtol=1e-3) +end + +n = 20 +x_sol = rand(1:Int(floor(n/4)), n) +diffi = x_sol + 0.3*rand([-1,1], n) + +@testset "Unit Simplex LMO" 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 + + N = sum(x_sol) + floor(n/2) + sblmo = Boscia.UnitSimplexSimpleBLMO(N) + + x, _, result = + Boscia.solve(f, grad!, sblmo, fill(0.0, n), fill(N, n), collect(1:n), n) + + @test sum(isapprox.(x, x_sol, atol=1e-6, rtol=1e-2)) == n + @test isapprox(f(x), f(result[:raw_solution]), atol=1e-6, rtol=1e-3) +end diff --git a/test/time_limit.jl b/test/time_limit.jl index 429ca7b96..65cafbb3a 100644 --- a/test/time_limit.jl +++ b/test/time_limit.jl @@ -63,5 +63,5 @@ time_limit = 30.0 time_taken = float(Dates.value(Dates.now() - start_time)) / 1000 @test sum(ai' * x) <= bi + 1e-6 @test f(x) <= f(result[:raw_solution]) + 1e-6 - @test time_taken <= time_limit + 10 + @test result[:total_time_in_sec] <= time_limit + 5 end