Indicator update best bound (#175)

* benchmark to compare v0.2.1 with v0.2.2

* update best bound v0.2.2

* using Cbc

* use Cbc at the right spot...

* branch variable dependent on objective

* have the change in `add_var!`

* changelog and v0.2.2

CHANGELOG.md

@@ -1,6 +1,10 @@
 # ConstrainSolver.jl - Changelog
 
-## v0.2.1 
+## v0.2.2 (26th of June 2020)
+- Actually use best bound [#175](https://github.com/Wikunia/ConstraintSolver.jl/pull/175)
+- Select next var based on objective (still hacky solution) [#176](https://github.com/Wikunia/ConstraintSolver.jl/issues/176)
+
+## v0.2.1 (26th of June 2020)
 - Bugfixes in indicator constraint [#170](https://github.com/Wikunia/ConstraintSolver.jl/issues/170)
   - Calling finished constraints and other functions for i.e `TableConstraint` as an inner constraint
   - Use correct best bound when inactive vs active

Project.toml

@@ -1,7 +1,7 @@
 name = "ConstraintSolver"
 uuid = "e0e52ebd-5523-408d-9ca3-7641f1cd1405"
 authors = ["Ole Kröger <[email protected]>"]
-version = "0.2.1"
+version = "0.2.2"
 
 [deps]
 Formatting = "59287772-0a20-5a39-b81b-1366585eb4c0"

benchmark/benchmarks.jl

@@ -1,6 +1,6 @@
 using BenchmarkTools
 using ConstraintSolver, JuMP, MathOptInterface
-using GLPK, JSON
+using GLPK, JSON, Cbc
 
 const CS = ConstraintSolver
 const MOI = MathOptInterface
@@ -23,6 +23,8 @@ SUITE["eternity"] = BenchmarkGroup(["alldifferent", "table", "equal"])
 # compiling run 
 solve_eternity("eternity_6x5"; height=6, width=5)
 SUITE["eternity"]["6x5"] = @benchmarkable solve_eternity("eternity_6x5"; height=6, width=5) seconds=30
+SUITE["eternity"]["5x5_opt"] = @benchmarkable solve_eternity("eternity_5x5"; height=5, width=5, optimize=true) seconds=120
+SUITE["eternity"]["5x5_opt_ind"] = @benchmarkable solve_eternity("eternity_5x5"; height=5, width=5, optimize=true, indicator=true) seconds=120
 SUITE["eternity"]["5x5_all"] = @benchmarkable solve_eternity("eternity_5x5"; all_solutions=true) seconds=30
 
 include(joinpath(dir, "benchmark/lp/benchmark.jl"))

benchmark/eternity/benchmark.jl

@@ -27,7 +27,7 @@ function get_rotations(puzzle)
     return rotations
 end
 
-function solve_eternity(fname="eternity_7"; height=nothing, width=nothing, all_solutions=false)
+function solve_eternity(fname="eternity_7"; height=nothing, width=nothing, all_solutions=false, optimize=false, indicator=false)
     puzzle = read_puzzle(fname)
     rotations = get_rotations(puzzle)
     npieces = size(puzzle)[1]
@@ -36,15 +36,27 @@ function solve_eternity(fname="eternity_7"; height=nothing, width=nothing, all_s
     ncolors = maximum(puzzle[:,2:end])
 
     m = Model(optimizer_with_attributes(CS.Optimizer, "logging" => [], "all_solutions"=>all_solutions))
+    if optimize
+        cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
+        m = Model(optimizer_with_attributes(CS.Optimizer, "logging" => [], "all_solutions"=>all_solutions, "lp_optimizer" => cbc_optimizer))
+    end
+
     @variable(m, 1 <= p[1:height, 1:width] <= npieces, Int)
     @variable(m, 0 <= pu[1:height, 1:width] <= ncolors, Int)
     @variable(m, 0 <= pr[1:height, 1:width] <= ncolors, Int)
     @variable(m, 0 <= pd[1:height, 1:width] <= ncolors, Int)
     @variable(m, 0 <= pl[1:height, 1:width] <= ncolors, Int)
+    if indicator 
+        @variable(m, b, Bin)
+    end
 
     @constraint(m, p[:] in CS.AllDifferentSet())
     for i=1:height, j=1:width
-        @constraint(m, [p[i,j], pu[i,j], pr[i,j], pd[i,j], pl[i,j]] in CS.TableSet(rotations))
+        if indicator
+            @constraint(m, b => {[p[i,j], pu[i,j], pr[i,j], pd[i,j], pl[i,j]] in CS.TableSet(rotations)})
+        else
+            @constraint(m, [p[i,j], pu[i,j], pr[i,j], pd[i,j], pl[i,j]] in CS.TableSet(rotations))
+        end
     end
 
     # borders
@@ -75,9 +87,19 @@ function solve_eternity(fname="eternity_7"; height=nothing, width=nothing, all_s
         @constraint(m, pr[i,j] == pl[i, j+1])
     end
 
-    if width == height
+    if !optimize && indicator
+        @constraint(m, b == 1)
+    end
+
+    if !optimize && width == height
         start_piece = findfirst(i->count(c->c == 0, puzzle[i,:]) == 2,1:npieces)
         @constraint(m, p[1,1] == start_piece)
+    elseif optimize
+        if indicator
+            @objective(m, Max, 1000*b + p[1,1] + p[1,2])
+        else
+            @objective(m, Max, p[1,1] + p[1,2])
+        end
     end
 
     optimize!(m)

benchmark/run_benchmarks.jl

@@ -33,7 +33,7 @@ end
 if isinteractive() == false
     args = parse_commandline()
     using PkgBenchmark
-    using ConstraintSolver
+    using ConstraintSolver, Cbc
     using GitHub, JSON, Statistics
 
     github_auth = GitHub.authenticate(ENV["GITHUB_AUTH"])

src/ConstraintSolver.jl

@@ -90,6 +90,7 @@ function add_var!(com::CS.CoM, from::Int, to::Int; fix = nothing)
     push!(com.search_space, var)
     push!(com.subscription, Int[])
     push!(com.bt_infeasible, 0)
+    push!(com.var_in_obj, false)
     return var
 end
 
@@ -488,13 +489,7 @@ function update_best_bound!(backtrack_obj::BacktrackObj, com::CS.CoM, constraint
     further_pruning = true
     feasible = true
     for constraint in constraints
-        relevant = false
-        for obj_index in com.objective.indices
-            if obj_index in constraint.std.indices
-                relevant = true
-                break
-            end
-        end
+        relevant = any(com.var_in_obj[i] for i in constraint.std.indices)
         if relevant
             feasible = prune_constraint!(
                 com,

src/MOI_wrapper/objective.jl

@@ -14,6 +14,7 @@ end
 
 function MOI.set(model::Optimizer, ::MOI.ObjectiveFunction, func::SVF)
     check_inbounds(model, func)
+    model.inner.var_in_obj[func.variable.value] = true
     model.inner.objective =
         SingleVariableObjective(func, func.variable.value, [func.variable.value])
     return
@@ -24,6 +25,7 @@ function MOI.set(model::Optimizer, ::MOI.ObjectiveFunction, func::SAF{T}) where
     indices = [func.terms[i].variable_index.value for i = 1:length(func.terms)]
     coeffs = [func.terms[i].coefficient for i = 1:length(func.terms)]
     lc = LinearCombination(indices, coeffs)
+    model.inner.var_in_obj[indices] .= true
     model.inner.objective = LinearCombinationObjective(func, lc, func.constant, indices)
     return
 end

src/MOI_wrapper/variables.jl

@@ -75,6 +75,7 @@ function MOI.add_variable(model::Optimizer)
     model.variable_info[index].changes = changes
     push!(model.inner.subscription, Int[])
     push!(model.inner.bt_infeasible, 0)
+    push!(model.inner.var_in_obj, false)
     addupd_var_in_inner_model(model, index)
     return MOI.VariableIndex(index)
 end

src/branching.jl

@@ -92,18 +92,29 @@ function get_next_branch_variable(com::CS.CoM)
     biggest_inf = -1
     best_ind = -1
     biggest_dependent = typemax(Int)
+    is_in_objective = false
     found = false
 
     for ind = 1:length(com.search_space)
         if !isfixed(com.search_space[ind])
             num_pvals = nvalues(com.search_space[ind])
             inf = com.bt_infeasible[ind]
-            if inf >= biggest_inf
-                if inf > biggest_inf || num_pvals < lowest_num_pvals
-                    lowest_num_pvals = num_pvals
-                    biggest_inf = inf
-                    best_ind = ind
-                    found = true
+            if !is_in_objective && com.var_in_obj[ind]
+                is_in_objective = true
+                lowest_num_pvals = num_pvals
+                biggest_inf = inf
+                best_ind = ind
+                found = true
+                continue
+            end
+            if !is_in_objective || com.var_in_obj[ind]
+                if inf >= biggest_inf
+                    if inf > biggest_inf || num_pvals < lowest_num_pvals
+                        lowest_num_pvals = num_pvals
+                        biggest_inf = inf
+                        best_ind = ind
+                        found = true
+                    end
                 end
             end
         end

src/constraints/indicator.jl

@@ -136,7 +136,6 @@ the possible values the table constraint allows. `var_idx`, `lb` and `ub` don't
 Additionally only a rough estimated bound is used which can be computed relatively fast. 
 This method calls the inner_constraint method if it exists and the indicator is activated.
 """
-#=
 function update_best_bound_constraint!(com::CS.CoM,
     constraint::IndicatorConstraint,
     fct::Union{MOI.VectorOfVariables, VAF{T}},
@@ -162,7 +161,6 @@ function update_best_bound_constraint!(com::CS.CoM,
     end
     return true
 end
-=#
 
 function single_reverse_pruning_constraint!(
     com::CoM,

src/type_inits.jl

@@ -88,7 +88,8 @@ function ConstraintSolverModel(::Type{T} = Float64) where {T<:Real}
         1, # c_backtrack_idx
         Vector{BacktrackObj{T}}(), # backtrack_vec
         MOI.FEASIBILITY_SENSE, #
-        NoObjective(), #
+        NoObjective(), # 
+        Vector{Bool}(), # var_in_obj
         get_traverse_strategy(), 
         get_branch_split(),
         zero(T), # best_sol,

src/types.jl

@@ -330,6 +330,7 @@ mutable struct ConstraintSolverModel{T<:Real}
     backtrack_vec::Vector{BacktrackObj{T}}
     sense::MOI.OptimizationSense
     objective::ObjectiveFunction
+    var_in_obj::Vector{Bool} # saves whether a variable is part of the objective function
     traverse_strategy::Val
     branch_split::Val
     best_sol::T # Objective of the best solution