Prototyping custom model/optimizer

Project.toml

@@ -10,6 +10,7 @@ DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 DynamicSparseArrays = "8086fd22-9a0c-46a5-a6c8-6e24676501fe"
+KnapsackLib = "86859df6-51c5-4863-9ac2-2c1ab8e53eb2"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
@@ -18,7 +19,6 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 
 [compat]
-BlockDecomposition = "~1.3"
 DataStructures = "0.17, 0.18"
 DynamicSparseArrays = "0.5"
 MathOptInterface = "~0.9.10"
@@ -30,7 +30,8 @@ julia = "1"
 ColunaDemos = "a54e61d4-7723-11e9-2469-af255fcaa246"
 GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
+KnapsackLib = "86859df6-51c5-4863-9ac2-2c1ab8e53eb2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["ColunaDemos", "GLPK", "JuMP", "Test"]
+test = ["ColunaDemos", "GLPK", "JuMP", "KnapsackLib", "Test"]

src/Algorithm/basic/solveipform.jl

@@ -77,14 +77,14 @@ end
 # SolveIpForm does not have child algorithms, therefore get_child_algorithms() is not defined
 
 # Dispatch on the type of the optimizer to return the parameters
-_optimizer_params(algo::SolveIpForm, ::MoiOptimizer) = algo.moi_params
-_optimizer_params(algo::SolveIpForm, ::UserOptimizer) = algo.user_params
-# TODO : custom optimizer
-_optimizer_params(::SolveIpForm, ::NoOptimizer) = nothing
+_optimizer_params(::Formulation, algo::SolveIpForm, ::MoiOptimizer) = algo.moi_params
+_optimizer_params(::Formulation, algo::SolveIpForm, ::UserOptimizer) = algo.user_params
+_optimizer_params(form::Formulation, algo::SolveIpForm, ::CustomOptimizer) = getinner(getoptimizer(form, algo.optimizer_id))
+_optimizer_params(::Formulation, ::SolveIpForm, ::NoOptimizer) = nothing
 
 function run!(algo::SolveIpForm, env::Env, form::Formulation, input::OptimizationInput)::OptimizationOutput
     opt = getoptimizer(form, algo.optimizer_id)
-    params = _optimizer_params(algo, opt)
+    params = _optimizer_params(form, algo, opt)
     if params !== nothing
         return run!(params, env, form, input; optimizer_id = algo.optimizer_id)
     end
@@ -99,7 +99,7 @@ run!(algo::SolveIpForm, env::Env, reform::Reformulation, input::OptimizationInpu
 ################################################################################
 function get_units_usage(algo::SolveIpForm, form::Formulation)
     opt = getoptimizer(form, algo.optimizer_id)
-    params = _optimizer_params(algo, opt)
+    params = _optimizer_params(form, algo, opt)
     if params !== nothing
         return get_units_usage(params, form)
     end
@@ -135,7 +135,8 @@ function get_units_usage(::UserOptimize, spform::Formulation{DwSp})
     return units_usage
 end
 
-# TODO : get_units_usage of CustomOptimize
+# no get_units_usage of CustomOptimize because it directly calls the
+# get_units_usage of the custom optimizer
 
 ################################################################################
 # run! methods (depends on the type of the optimizer)
@@ -270,4 +271,5 @@ function run!(
     return OptimizationOutput(result)
 end
 
-# TODO : run! of CustomOptimize
+# No run! method for CustomOptimize because it directly calls the run! method
+# of the custom optimizer

src/MathProg/MathProg.jl

@@ -43,10 +43,10 @@ include("MOIinterface.jl")
 
 # TODO : clean up
 # Types
-export  MaxSense, MinSense, MoiOptimizer,
+export  MaxSense, MinSense,
         Id, ConstrSense, VarSense,
         FormId, FormulationPhase, Annotations,
-        Counter, UserOptimizer, NoOptimizer, MoiObjective
+        Counter, MoiObjective
 
 # Methods
 export no_optimizer_builder, set_original_formulation!,
@@ -110,4 +110,7 @@ export PrimalBound, DualBound, PrimalSolution, DualSolution, ObjValues,
 # Methods related to projections
 export projection_is_possible, proj_cols_on_rep
 
+# Optimizers of formulations
+export MoiOptimizer, CustomOptimizer, UserOptimizer, NoOptimizer
+
 end

src/MathProg/optimizerwrappers.jl

@@ -147,3 +147,12 @@ function write_to_LP_file(form::Formulation, optimizer::MoiOptimizer, filename::
     MOI.copy_to(dest, src)
     MOI.write_to_file(dest, filename)
 end
+
+"""
+    CustomOptimizer <: AbstractOptimizer
+"""
+struct CustomOptimizer <: AbstractOptimizer
+    inner::BD.AbstractCustomOptimizer
+end
+
+getinner(optimizer::CustomOptimizer) = optimizer.inner

src/decomposition.jl

@@ -427,6 +427,7 @@ end
 
 _optimizerbuilder(opt::Function) = () -> UserOptimizer(opt)
 _optimizerbuilder(opt::MOI.AbstractOptimizer) = () -> MoiOptimizer(opt)
+_optimizerbuilder(opt::BD.AbstractCustomOptimizer) = () -> CustomOptimizer(opt)
 
 function getoptimizerbuilders(prob::Problem, ann::BD.Annotation)
     optimizers = BD.getoptimizerbuilders(ann)

test/Project.toml

@@ -8,4 +8,5 @@ MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 LightGraphs = "093fc24a-ae57-5d10-9952-331d41423f4d"
-Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
+Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
+KnapsackLib = "86859df6-51c5-4863-9ac2-2c1ab8e53eb2"

test/interfaces/model.jl

@@ -0,0 +1,117 @@
+# In this test, we use the Martinelli's knapsack solver pkg ( https://github.com/rafaelmartinelli/KnapsackLib.jl)
+# to test the interface of custom models/solvers.
+
+using KnapsackLib
+mutable struct KnapsackLibModel <: Coluna.MathProg.AbstractFormulation
+    nbitems::Int
+    costs::Vector{Float64}
+    weights::Vector{Float64}
+    capacity::Float64
+    job_to_jumpvar::Dict{Int, JuMP.VariableRef}
+    #varids::Vector{Coluna.MathProg.VarId}
+    #map::Dict{Coluna.MathProg.VarId,Float64}
+end
+KnapsackLibModel(nbitems) = KnapsackLibModel(
+    nbitems, zeros(Float64, nbitems), zeros(Float64, nbitems), 0.0,
+    Dict{Int, JuMP.VariableRef}()
+)
+setcapacity!(model::KnapsackLibModel, cap) = model.capacity = cap
+setweight!(model::KnapsackLibModel, j::Int, w) = model.weights[j] = w
+setcost!(model::KnapsackLibModel, j::Int, c) = model.costs[j] = c
+map!(model::KnapsackLibModel, j::Int, x::JuMP.VariableRef) = model.job_to_jumpvar[j] = x
+
+coluna_backend(model::MOI.Utilities.CachingOptimizer) = coluna_backend(model.optimizer)
+coluna_backend(b::MOI.Bridges.AbstractBridgeOptimizer) = coluna_backend(b.model)
+coluna_backend(model) = model
+
+function get_coluna_varid(model::KnapsackLibModel, form, j::Int)
+    jumpvar = model.job_to_jumpvar[j]
+    opt = coluna_backend(backend(jumpvar.model))
+    return Coluna._get_orig_varid_in_form(opt, form, jumpvar.index)
+end
+mutable struct KnapsackLibOptimizer <: BlockDecomposition.AbstractCustomOptimizer
+    model::KnapsackLibModel
+end
+
+function Coluna.Algorithm.get_units_usage(opt::KnapsackLibOptimizer, form) # form is Coluna Formulation
+    println("\e[41m get units usage \e[00m")
+    units_usage = Tuple{AbstractModel, Coluna.ColunaBase.UnitType, Coluna.ColunaBase.UnitAccessMode}[]
+    # TODO : the abstract model is KnapsackLibModel (opt.model)
+    return units_usage
+end
+
+function _rfl(val::Float64)::Integer
+    rf_val = Integer(floor(val + val * 1e-10 + 1e-6))
+    rf_val += rf_val < val - 1 + 1e-6 ? 1 : 0
+    return rf_val
+end
+
+function _scale_to_int(vals...)
+    max_val = maximum(vals)
+    scaling_factor = typemax(Int) / (length(vals) + 1) / max_val
+    return map(x -> _rfl(scaling_factor * x), vals)
+end
+
+function Coluna.Algorithm.run!(
+    opt::KnapsackLibOptimizer, ::Coluna.Env, form::Coluna.MathProg.Formulation,
+    input::Coluna.Algorithm.OptimizationInput; kw...
+)
+    ws = _scale_to_int(opt.model.capacity, opt.model.weights...)
+    cs = _scale_to_int(opt.model.costs...)
+    items = [KnapItem(w,c) for (w,c) in zip(ws[2:end], cs)]
+    data = KnapData(ws[1], items)
+    _, selected = solveKnapExpCore(data)
+    optimal = sum(opt.model.costs[j] for j in selected)
+    @show optimal
+    @show selected
+    for j in selected
+        @show Coluna.MathProg.getname(form, get_coluna_varid(opt.model, form, j))
+    end
+
+    error("run! method of custom optimizer reached !")
+    return
+end
+
+################################################################################
+# User model
+################################################################################
+function knpcustommodel()
+    @testset "knapsack custom model" begin
+        data = CLD.GeneralizedAssignment.data("play2.txt")
+        coluna = JuMP.optimizer_with_attributes(
+            Coluna.Optimizer,
+            "params" => CL.Params(solver = ClA.TreeSearchAlgorithm()),
+            "default_optimizer" => GLPK.Optimizer
+        )
+
+        model = BlockModel(coluna; direct_model = true)
+        @axis(M, data.machines)
+        @variable(model, x[m in M, j in data.jobs], Bin)
+        @constraint(model, 
+            sp[j in data.jobs], sum(x[m,j] for m in data.machines) == 1
+        )
+        @objective(model, Min,
+            sum(data.cost[j,m]*x[m,j] for m in M, j in data.jobs)
+        )
+
+        @dantzig_wolfe_decomposition(model, dec, M)
+
+        sp = getsubproblems(dec)
+        for m in M
+            knp_model = KnapsackLibModel(length(data.jobs))
+            setcapacity!(knp_model, data.capacity[m])
+            for j in data.jobs
+                setweight!(knp_model, j, data.weight[j,m])
+                setcost!(knp_model, j, data.cost[j,m])
+                map!(knp_model, j, x[m,j])
+            end
+            knp_optimizer = KnapsackLibOptimizer(knp_model)
+            specify!(sp[m], solver = knp_optimizer) ##model = knp_model)
+        end
+
+        optimize!(model)
+    end
+    exit()
+end
+
+knpcustommodel()

test/runtests.jl

@@ -23,6 +23,7 @@ const ClA = Coluna.Algorithm
 
 include("unit/unit_tests.jl")
 include("MathOptInterface/MOI_wrapper.jl")
+include("interfaces/model.jl")
 include("issues_tests.jl")
 include("show_functions_tests.jl")
 include("full_instances_tests.jl")