Merge pull request #269 from scipopt/heuristic

Heuristic plugin

.codecov.yml

@@ -2,3 +2,4 @@ comment: false
 ignore:
   - "**/compat.jl"
   - "**/wrapper/*.jl"
+  - "**LibSCIP.jl"

src/MOI_wrapper.jl

@@ -57,6 +57,7 @@ mutable struct Optimizer <: MOI.AbstractOptimizer
             Dict(),
             Dict(),
             Dict(),
+            Dict(),
             [],
         )
 
@@ -238,7 +239,7 @@ function MOI.empty!(o::Optimizer)
     @SCIP_CALL SCIP.SCIPcreateProbBasic(scip[], "")
     # create a new problem
     o.inner =
-        SCIPData(scip, Dict(), Dict(), 0, 0, Dict(), Dict(), Dict(), Dict(), Dict(), [])
+        SCIPData(scip, Dict(), Dict(), 0, 0, Dict(), Dict(), Dict(), Dict(), Dict(), Dict(), [])
     # reapply parameters
     for pair in o.params
         set_parameter(o.inner, pair.first, pair.second)
@@ -371,3 +372,4 @@ include(joinpath("MOI_wrapper", "objective.jl"))
 include(joinpath("MOI_wrapper", "results.jl"))
 include(joinpath("MOI_wrapper", "conshdlr.jl"))
 include(joinpath("MOI_wrapper", "sepa.jl"))
+include(joinpath("MOI_wrapper", "heuristic.jl"))

src/MOI_wrapper/heuristic.jl

@@ -0,0 +1,29 @@
+
+function include_heuristic(
+    o::Optimizer,
+    heuristic::HT;
+    name="",
+    description="",
+    dispchar='_',
+    priority=10000,
+    frequency=1,
+    frequency_offset=0,
+    maximum_depth=-1,
+    timing_mask=SCIP_HEURTIMING_BEFORENODE,
+    usessubscip=false,
+) where {HT}
+    return include_heuristic(
+        o.inner.scip[],
+        heuristic,
+        o.inner.heuristic_storage;
+        name=name,
+        description=description,
+        dispchar=dispchar,
+        priority=priority,
+        frequency=frequency,
+        frequency_offset=frequency_offset,
+        maximum_depth=maximum_depth,
+        timing_mask=timing_mask,
+        usessubscip=usessubscip,
+    )
+end

src/SCIP.jl

@@ -3,6 +3,7 @@ module SCIP
 # assorted utility functions
 include("util.jl")
 
+
 # load deps, version check
 include("init.jl")
 
@@ -12,13 +13,12 @@ include("wrapper.jl")
 # memory management
 include("scip_data.jl")
 
-# separators
 include("sepa.jl")
 
-# cut selectors
 include("cut_selector.jl")
 
-# branching rule
+include("heuristic.jl")
+
 include("branching_rule.jl")
 
 # constraint handlers

src/heuristic.jl

@@ -0,0 +1,169 @@
+# heuristic interface
+# it is recommended to check https://scipopt.org/doc/html/HEUR.php for key concepts and interface
+
+"""
+Abstract class for Heuristic.
+A heuristic must implement `find_primal_solution`.
+It also stores all user data that must be available to run the heuristic.
+"""
+abstract type Heuristic end
+
+"""
+    find_primal_solution(
+        scip::Ptr{SCIP_},
+        heur::Heuristic,
+        heurtiming::Heurtiming,
+        nodeinfeasible::Bool,
+        heur_ptr::Ptr{SCIP_HEUR},
+    ) -> (retcode, result, solutions)
+
+It must attempt to find primal solution(s).
+`retcode` indicates whether the selection went well.
+A typical result would be `SCIP_SUCCESS`, and retcode `SCIP_OKAY`.
+`solutions` is a vector of added SCIP_SOL pointers.
+Use the methods `create_scipsol` and `SCIPsetSolVal` to build solutions.
+Do not add them to SCIP directly (i.e. do not call `SCIPtrySolFree`).
+"""
+function find_primal_solution(scip, heur, heurtiming, nodeinfeasible, heur_ptr) end
+
+function _find_primal_solution_callback(
+    scip::Ptr{SCIP_},
+    heur_::Ptr{SCIP_HEUR},
+    heurtiming::SCIP_HEURTIMING,
+    nodeinfeasible_::SCIP_Bool,
+    result_::Ptr{SCIP_RESULT},
+)
+    heurdata::Ptr{SCIP_HEURDATA} = SCIPheurGetData(heur_)
+    heur = unsafe_pointer_to_objref(heurdata)
+    nodeinfeasible = nodeinfeasible_ == SCIP.TRUE
+    (retcode, result, solutions) = find_primal_solution(
+        scip,
+        heur,
+        heurtiming,
+        nodeinfeasible,
+        heur_,
+    )::Tuple{SCIP_RETCODE,SCIP_RESULT,Vector{Ptr{SCIP_SOL}}}
+    if retcode != SCIP_OKAY
+        return retcode
+    end
+    if result == SCIP_FOUNDSOL
+        @assert length(solutions) > 0
+    end
+    found_solution = false
+    for sol in solutions
+        stored = Ref{SCIP_Bool}(SCIP.FALSE)
+        @SCIP_CALL SCIPtrySolFree(
+            scip,
+            Ref(sol),
+            SCIP.FALSE,
+            SCIP.FALSE,
+            SCIP.TRUE,
+            SCIP.TRUE,
+            SCIP.TRUE,
+            stored, 
+	    )
+        if stored[] != SCIP.TRUE
+            @warn "Primal solution not feasible"
+        else
+            found_solution = true
+        end
+    end
+    result = if found_solution
+        SCIP_FOUNDSOL
+    else
+        SCIP_DIDNOTFIND
+    end
+
+    unsafe_store!(result_, result)
+    return retcode
+end
+
+function _heurfree(::Ptr{SCIP_}, heur::Ptr{SCIP_HEUR})
+    # just like sepa, free the data on the SCIP side,
+    # the Julia GC will take care of the objects
+    SCIPheurSetData(heur, C_NULL)
+    return SCIP_OKAY
+end
+
+"""
+Includes a heuristic plugin in SCIP and stores it in heuristic_storage.
+"""
+function include_heuristic(
+    scip::Ptr{SCIP_},
+    heuristic::HT,
+    heuristic_storage::Dict{Any,Ptr{SCIP_HEUR}};
+    name="",
+    description="",
+    dispchar='_',
+    priority=10000,
+    frequency=1,
+    frequency_offset=0,
+    maximum_depth=-1,
+    timing_mask=SCIP_HEURTIMING_BEFORENODE,
+    usessubscip=false,
+) where {HT<:Heuristic}
+
+    # ensure a unique name for the cut selector
+    if name == ""
+        name = "heuristic_$(string(HT))"
+    end
+    if dispchar == '_'
+        dispchar = name[end]
+    end
+
+    heur__ = Ref{Ptr{SCIP_HEUR}}(C_NULL)
+    if !ismutable(heuristic)
+        throw(
+            ArgumentError("The heuristic structure must be a mutable type"),
+        )
+    end
+
+    heurdata_ = pointer_from_objref(heuristic)
+    heur_callback = @cfunction(
+        _find_primal_solution_callback,
+        SCIP_RETCODE,
+        (
+            Ptr{SCIP_},
+            Ptr{SCIP_HEUR},
+            SCIP_HEURTIMING,
+            SCIP_Bool,
+            Ptr{SCIP_RESULT},
+        ),
+    )
+    @SCIP_CALL SCIPincludeHeurBasic(
+        scip,
+        heur__,
+        name,
+        description,
+        dispchar,
+        priority,
+        frequency,
+        frequency_offset,
+        maximum_depth,
+        timing_mask,
+        usessubscip,
+        heur_callback,
+        heurdata_,
+    )
+
+    @assert heur__[] != C_NULL
+
+    @SCIP_CALL SCIPsetHeurFree(
+        scip,
+        heur__[],
+        @cfunction(_heurfree, SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_HEUR})),
+    )
+
+    # store heuristic in storage (avoids GC-ing it)
+    heuristic_storage[heuristic] = heur__[]
+end
+
+"""
+create_scipsol(scip::Ptr{SCIP_}, heur_::Ptr{SCIP_HEUR}) -> Ptr{SCIP_SOL}
+Convenience wrapper to create a 
+"""
+function create_scipsol(scip::Ptr{SCIP_}, heur_::Ptr{SCIP_HEUR})
+    sol__ = Ref{Ptr{SCIP_SOL}}(C_NULL)
+    @SCIP_CALL SCIPcreateSol(scip, sol__, heur_)
+    return sol__[]
+end

src/scip_data.jl

@@ -46,7 +46,8 @@ mutable struct SCIPData
     # User-defined cut selectors and branching rules
     cutsel_storage::Dict{Any,Ptr{SCIP_CUTSEL}}
     branchrule_storage::Dict{Any,Ptr{SCIP_BRANCHRULE}}
-
+    heuristic_storage::Dict{Any,Ptr{SCIP_HEUR}}
+
     # to store expressions for release
     nonlinear_storage::Vector{NonlinExpr}
 end

test/cutcallback.jl

@@ -100,10 +100,7 @@ end
     MOI.set(
         optimizer,
         MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
-        MOI.ScalarAffineFunction(
-            MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]),
-            0.0,
-        ),
+        1.0 * x + 1.0 * y,
     )
     MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)
 
@@ -112,7 +109,7 @@ end
         MOI.submit(
             optimizer,
             MOI.UserCut{SCIP.CutCbData}(cb_data),
-            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0], [x]), 0.0),
+            1.0 * x,
             MOI.LessThan(0.0),
         )
         calls += 1

test/heuristic.jl

@@ -0,0 +1,63 @@
+import MathOptInterface as MOI
+using SCIP
+using LinearAlgebra
+using Test
+
+mutable struct ZeroHeuristic <: SCIP.Heuristic
+end
+
+function SCIP.find_primal_solution(scip, ::ZeroHeuristic, heurtiming, nodeinfeasible::Bool, heur_ptr)
+    @assert SCIP.SCIPhasCurrentNodeLP(scip) == SCIP.TRUE
+    result = SCIP.SCIP_DIDNOTRUN
+    sol = SCIP.create_scipsol(scip, heur_ptr)
+    vars = SCIP.SCIPgetVars(scip)
+    nvars = SCIP.SCIPgetNVars(scip)
+    var_vec = unsafe_wrap(Array, vars, nvars)
+    for var in var_vec
+        SCIP.@SCIP_CALL SCIP.SCIPsetSolVal(
+            scip,
+            sol,
+            var,
+            0.0, 
+        )
+    end
+    result = SCIP.SCIP_SUCCESS
+    return (SCIP.SCIP_OKAY, result, [sol])
+end
+
+@testset "Basic heuristic properties" begin
+    o = SCIP.Optimizer(; presolving_maxrounds=0)
+    name = "zero_heuristic"
+    description = "description"
+    priority = 1
+    heur = ZeroHeuristic()
+    SCIP.include_heuristic(o, heur, name=name, description=description, priority=priority)
+
+    heur_pointer = o.inner.heuristic_storage[heur]
+    @test unsafe_string(SCIP.LibSCIP.SCIPheurGetName(heur_pointer)) == name
+    @test unsafe_string(SCIP.LibSCIP.SCIPheurGetDesc(heur_pointer)) ==
+          description
+    @test SCIP.LibSCIP.SCIPheurGetPriority(heur_pointer) == priority
+    @test SCIP.LibSCIP.SCIPheurGetData(heur_pointer) ==
+          pointer_from_objref(heur)
+
+    x = MOI.add_variables(o, 10)
+    MOI.add_constraint.(o, x, MOI.Integer())
+    MOI.add_constraint.(o, x, MOI.GreaterThan(-0.1))
+    MOI.add_constraint.(o, x, MOI.LessThan(2.3))
+    MOI.add_constraint(o, sum(x; init=0.0), MOI.LessThan(12.5))
+    for _ in 1:5
+        MOI.add_constraint(
+            o,
+            2.0 * dot(rand(10), x),
+            MOI.LessThan(10.0 + 2 * rand()),
+        )
+    end
+    func = -dot(rand(10), x)
+    MOI.set(o, MOI.ObjectiveFunction{typeof(func)}(), func)
+    MOI.set(o, MOI.ObjectiveSense(), MOI.MIN_SENSE)
+
+    MOI.optimize!(o)
+    @test MOI.get(o, MOI.TerminationStatus()) == MOI.OPTIMAL
+
+end

test/runtests.jl

@@ -44,6 +44,9 @@ end
 @testset "branching rule" begin
     include("branchrule.jl")
 end
+@testset "heuristic" begin
+    include("heuristic.jl")
+end
 
 const MOI_BASE_EXCLUDED = [
     "Indicator_LessThan", # indicator must be binary error in SCIP