Implementing #57: Real coefficients (#58)

* Implementing #57: Real coefficients for eq_sum constraint
* `SolverOptions()` constructor
* `LinearVariables` -> `LinearCombination`

Author: Wikunia

benchmark/killer_sudoku/cs.jl

@@ -1,4 +1,4 @@
-using ConstraintSolver, MathOptInterface, JSON
+using ConstraintSolver, JuMP, MathOptInterface, JSON
 
 if !@isdefined CS 
     const CS = ConstraintSolver
@@ -57,8 +57,8 @@ function solve_all(filenames; benchmark=false, single_times=true)
             GC.enable(true)
         end
         if !benchmark
-            @show m.inner.info
             println("Status: ", status)
+            @show m.inner.info
             solution = zeros(Int, 9, 9)
             for r=1:9
                 solution[r,:] = [MOI.get(m, MOI.VariablePrimal(), x[r][c][1]) for c=1:9]

benchmark/sudoku/cs.jl

@@ -1,6 +1,8 @@
-using ConstraintSolver, MathOptInterface
+using ConstraintSolver, JuMP, MathOptInterface
 
-const CS = ConstraintSolver
+if !@isdefined CS 
+    const CS = ConstraintSolver
+end
 const MOI = MathOptInterface
 const MOIU = MOI.Utilities
 include("../../test/sudoku_fcts.jl")
@@ -57,8 +59,8 @@ function solve_all(grids; benchmark=false, single_times=true)
             GC.enable(true)
         end
         if !benchmark
-            @show m.inner.info
             println("Status: ", status)
+            @show m.inner.info
             solution = zeros(Int, 9, 9)
             for r=1:9
                 solution[r,:] = [MOI.get(m, MOI.VariablePrimal(), x[r][c][1]) for c=1:9]

src/ConstraintSolver.jl

@@ -83,26 +83,60 @@ struct NotEqualSet{T} <: MOI.AbstractScalarSet
     value :: T
 end
 
-mutable struct LinearVariables
+mutable struct LinearCombination{T <: Real}
     indices             :: Vector{Int}
-    coeffs              :: Vector{Int}
+    coeffs              :: Vector{T}
 end
 
-mutable struct LinearConstraint <: Constraint
+mutable struct LinearConstraint{T <: Real} <: Constraint
     idx                 :: Int
     fct                 :: Function
     indices             :: Vector{Int}
     pvals               :: Vector{Int}
-    coeffs              :: Vector{Int}
+    coeffs              :: Vector{T}
     operator            :: Symbol
-    rhs                 :: Int
+    rhs                 :: T
     in_all_different    :: Bool
-    mins                :: Vector{Int}
-    maxs                :: Vector{Int}
-    pre_mins            :: Vector{Int}
-    pre_maxs            :: Vector{Int}
+    mins                :: Vector{T}
+    maxs                :: Vector{T}
+    pre_mins            :: Vector{T}
+    pre_maxs            :: Vector{T}
     hash                :: UInt64
-    LinearConstraint() = new()
+end
+
+function LinearConstraint(fct::Function, operator::Symbol, indices::Vector{Int}, coeffs::Vector{T}, rhs::Real) where T <: Real
+    # get common type for rhs and coeffs
+    promote_T = promote_type(typeof(rhs), eltype(coeffs))
+    if promote_T != eltype(coeffs)
+        coeffs = convert.(promote_T, coeffs)
+    end
+    if promote_T != typeof(rhs)
+        rhs = convert(promote_T, rhs)
+    end
+    maxs = zeros(promote_T, length(indices))
+    mins = zeros(promote_T, length(indices))
+    pre_maxs = zeros(promote_T, length(indices))
+    pre_mins = zeros(promote_T, length(indices))
+    # this can be changed later in `set_in_all_different!` but needs to be initialized with false
+    in_all_different = false
+    pvals = Int[]
+
+    lc = LinearConstraint(
+        0, # idx will be filled later
+        fct,
+        indices,
+        pvals, 
+        coeffs,
+        operator,
+        rhs,
+        in_all_different,
+        mins,
+        maxs,
+        pre_mins,
+        pre_maxs,
+        zero(UInt64)
+    )
+    return lc
 end
 
 mutable struct BacktrackObj
@@ -147,10 +181,12 @@ mutable struct ConstraintSolverModel
     info                :: CSInfo
     input               :: Dict{Symbol,Any}
     logs                :: Vector{TreeLogNode}
+    options             :: SolverOptions
 end
 
 const CoM = ConstraintSolverModel
 
+include("util.jl")
 include("MOI_wrapper/MOI_wrapper.jl")
 include("printing.jl")
 include("logs.jl")
@@ -185,7 +221,8 @@ function ConstraintSolverModel()
         Vector{Int}(), # solutions
         CSInfo(0, false, 0, 0, 0), # info
         Dict{Symbol,Any}(), # input
-        Vector{TreeLogNode}() # logs
+        Vector{TreeLogNode}(), # logs
+        SolverOptions() # options
     )
 end
 
@@ -629,9 +666,38 @@ function update_best_bound!(backtrack_obj::BacktrackObj, com::CS.CoM, constraint
     if backtrack_obj.best_bound != new_bb
         further_pruning = false
     end
+    if backtrack_obj.best_bound == com.best_bound
+        backtrack_obj.best_bound = new_bb
+        update_best_bound!(com)
+    else 
+        backtrack_obj.best_bound = new_bb
+    end
     return true, further_pruning
 end
 
+function update_best_bound!(com::CS.CoM)
+    if com.sense == MOI.MIN_SENSE
+        max_val = typemax(Int64)
+        com.best_bound = minimum([bo.status == :Open ? bo.best_bound : max_val for bo in com.backtrack_vec])
+    elseif com.sense == MOI.MAX_SENSE
+        min_val = typemin(Int64)
+        com.best_bound = maximum([bo.status == :Open ? bo.best_bound : min_val for bo in com.backtrack_vec])
+    else
+        com.best_bound = 0
+    end 
+end
+
+function set_state_to_best_sol!(com::CS.CoM, last_backtrack_id::Int)
+    obj_factor = com.sense == MOI.MIN_SENSE ? 1 : -1
+    backtrack_vec = com.backtrack_vec
+    # find one of the best solutions
+    sol, sol_id = findmin([backtrack_vec[sol_id].best_bound*obj_factor for sol_id in com.solutions])
+    backtrack_id = com.solutions[sol_id]
+    checkout_from_to!(com, last_backtrack_id, backtrack_id)
+    # prune the last step as checkout_from_to! excludes the to part
+    prune!(com, [backtrack_id])
+end
+
 """
     backtrack!(com::CS.CoM, max_bt_steps; sorting=true)
 
@@ -731,15 +797,7 @@ function backtrack!(com::CS.CoM, max_bt_steps; sorting=true)
         if l <= 0 || l > length(backtrack_vec)
             break
         end
-        if com.sense == MOI.MIN_SENSE
-            max_val = typemax(Int64)
-            com.best_bound = minimum([bo.status == :Open ? bo.best_bound : max_val for bo in backtrack_vec])
-        elseif com.sense == MOI.MAX_SENSE
-            min_val = typemin(Int64)
-            com.best_bound = maximum([bo.status == :Open ? bo.best_bound : min_val for bo in backtrack_vec])
-        else
-            com.best_bound = 0
-        end
+        update_best_bound!(com)
 
         ind = backtrack_obj.variable_idx
 
@@ -782,7 +840,6 @@ function backtrack!(com::CS.CoM, max_bt_steps; sorting=true)
         # first update the best bound (only constraints which have an index in the objective function)
         if com.sense != MOI.FEASIBILITY_SENSE
             feasible, further_pruning = update_best_bound!(backtrack_obj, com, constraints)
-
             if !feasible
                 com.info.backtrack_reverses += 1
                 continue
@@ -792,7 +849,7 @@ function backtrack!(com::CS.CoM, max_bt_steps; sorting=true)
         if further_pruning
             # prune completely start with all that changed by the fix or by updating best bound
             feasible = prune!(com)
-
+        
             if !feasible
                 com.info.backtrack_reverses += 1
                 continue
@@ -804,35 +861,36 @@ function backtrack!(com::CS.CoM, max_bt_steps; sorting=true)
         # no index found => solution found
         if !found
             new_sol = get_best_bound(com)
-            if length(com.solutions) == 0
+            if length(com.solutions) == 0 || obj_factor*new_sol <= obj_factor*com.best_sol
+                push!(com.solutions, backtrack_obj.idx)
                 com.best_sol = new_sol
-            elseif obj_factor*new_sol < obj_factor*com.best_sol
-                com.best_sol = new_sol
-            end
-            push!(com.solutions, backtrack_obj.idx)
-            if com.best_sol == com.best_bound
-                return :Solved
-            else
+                if com.best_sol == com.best_bound
+                    return :Solved
+                end
                 # set all nodes to :Worse if they can't achieve a better solution
                 for bo in backtrack_vec
                     if bo.status == :Open && obj_factor*bo.best_bound >= com.best_sol
                         bo.status = :Worse
                     end
                 end
                 continue
+            else
+                if com.best_sol == com.best_bound
+                    set_state_to_best_sol!(com, last_backtrack_id)
+                    return :Solved
+                end
+                continue
             end
         end
 
         if com.info.backtrack_fixes > max_bt_steps
             return :NotSolved
         end
 
-
-
         if com.input[:logs]
             com.logs[backtrack_obj.idx] = log_one_node(com, length(com.search_space), backtrack_obj.idx, step_nr)
         end
-
+        
         pvals = reverse!(values(com.search_space[ind]))
         last_backtrack_obj = backtrack_vec[last_backtrack_id]
         for pval in pvals
@@ -861,12 +919,7 @@ function backtrack!(com::CS.CoM, max_bt_steps; sorting=true)
         end
     end
     if length(com.solutions) > 0
-        # find one of the best solutions
-        sol, sol_id = findmin([backtrack_vec[sol_id].best_bound*obj_factor for sol_id in com.solutions])
-        backtrack_id = com.solutions[sol_id]
-        checkout_from_to!(com, last_backtrack_id, backtrack_id)
-        # prune the last step as checkout_from_to! excludes the to part
-        prune!(com, [backtrack_id])
+        set_state_to_best_sol!(com, last_backtrack_id)
         return :Solved
     end
     return :Infeasible
@@ -997,6 +1050,7 @@ function solve!(com::CS.CoM, options::SolverOptions)
         return :Solved
     end
 
+    com.options = options
     backtrack = options.backtrack
     max_bt_steps = options.max_bt_steps
     backtrack_sorting = options.backtrack_sorting

src/MOI_wrapper/constraints.jl

@@ -23,21 +23,14 @@ function MOI.add_constraint(model::Optimizer, func::SAF, set::MOI.EqualTo{Float6
         fix!(model.inner, model.variable_info[func.terms[1].variable_index.value], convert(Int64, set.value/func.terms[1].coefficient))
         return MOI.ConstraintIndex{SAF, MOI.EqualTo{Float64}}(0)
     end
-
-    lc = LinearConstraint()
+   
     indices = [v.variable_index.value for v in func.terms]
     coeffs = [v.coefficient for v in func.terms]
-    lc.fct = eq_sum
-    lc.indices = indices
-    lc.coeffs = coeffs
-    lc.operator = :(==)
-    lc.rhs = set.value
-    lc.maxs = zeros(Int, length(indices))
-    lc.mins = zeros(Int, length(indices))
-    lc.pre_maxs = zeros(Int, length(indices))
-    lc.pre_mins = zeros(Int, length(indices))
-    # this can be changed later in `set_in_all_different!` but needs to be initialized with false
-    lc.in_all_different = false
+    fct = eq_sum
+    operator = :(==)
+    rhs = set.value
+    
+    lc = LinearConstraint(fct, operator, indices, coeffs, rhs)
     lc.idx = length(model.inner.constraints)+1
 
     push!(model.inner.constraints, lc)