Retrieve disaggregated solution

src/MOIwrapper.jl

@@ -29,6 +29,7 @@ mutable struct Optimizer <: MOI.AbstractOptimizer
     constrs_on_single_var_to_names::Dict{MOI.ConstraintIndex, String}
     names_to_constrs::Dict{String, MOI.ConstraintIndex}
     result::OptimizationState
+    disagg_result::Union{Nothing, OptimizationState}
     default_optimizer_builder::Union{Nothing, Function}
 
     feasibility_sense::Bool # Coluna supports only Max or Min.
@@ -47,6 +48,7 @@ mutable struct Optimizer <: MOI.AbstractOptimizer
         model.constrs_on_single_var_to_names = Dict{MOI.ConstraintIndex, String}()
         model.names_to_constrs = Dict{String, MOI.ConstraintIndex}()
         model.result = OptimizationState(get_optimization_target(model.inner))
+        model.disagg_result = nothing
         model.default_optimizer_builder = nothing
         model.feasibility_sense = false
         return model
@@ -95,7 +97,7 @@ end
 MOI.get(optimizer::Coluna.Optimizer, ::MOI.SolverName) = "Coluna"
 
 function MOI.optimize!(optimizer::Optimizer)
-    optimizer.result = optimize!(
+    optimizer.result, optimizer.disagg_result = optimize!(
         optimizer.env, optimizer.inner, optimizer.annotations
     )
     return
@@ -588,9 +590,36 @@ function MOI.empty!(model::Coluna.Optimizer)
         set_default_optimizer_builder!(model.inner, model.default_optimizer_builder)
     end
     model.result = OptimizationState(get_optimization_target(model.inner))
+    model.disagg_result = nothing
     return
 end
 
+mutable struct ColumnInfo <: BD.AbstractColumnInfo
+    optimizer::Coluna.Optimizer
+    column_var_id::VarId
+    column_val::Float64
+end
+
+function BD.getsolutions(model::Coluna.Optimizer, k)
+    ip_primal_sol = model.disagg_result.ip_primal_sols[k]
+    sp_columns_info = Vector{ColumnInfo}()
+    for (varid, val) in ip_primal_sol
+        if getduty(varid) <= MasterCol
+            push!(sp_columns_info, ColumnInfo(model, varid, val))
+        end
+    end
+    return sp_columns_info
+end
+
+BD.value(info::ColumnInfo) = info.column_val
+
+function BD.value(info::ColumnInfo, index::MOI.VariableIndex)
+    varid = info.optimizer.env.varids[index]
+    origin_form_uid = getoriginformuid(info.column_var_id)
+    spform = get_dw_pricing_sps(info.optimizer.inner.re_formulation)[origin_form_uid]
+    return info.column_val * getprimalsolmatrix(spform)[varid, info.column_var_id]
+end
+
 function MOI.get(model::Coluna.Optimizer, ::MOI.NumberOfVariables)
     orig_form = get_original_formulation(model.inner)
     return length(getvars(orig_form))

src/optimize.jl

@@ -60,7 +60,7 @@ function optimize!(env::Env, prob::MathProg.Problem, annotations::Annotations)
     @logmsg LogLevel(-1) env.params
 
     TO.@timeit _to "Coluna" begin
-        optstate = optimize!(get_optimization_target(prob), env, init_pb, init_db)
+        outstate, algstate = optimize!(get_optimization_target(prob), env, init_pb, init_db)
     end
 
     env.kpis.elapsed_optimization_time = elapsed_optim_time(env)
@@ -71,9 +71,9 @@ function optimize!(env::Env, prob::MathProg.Problem, annotations::Annotations)
     TO.reset_timer!(_to)
 
     @logmsg LogLevel(0) "Terminated"
-    @logmsg LogLevel(0) string("Primal bound: ", get_ip_primal_bound(optstate))
-    @logmsg LogLevel(0) string("Dual bound: ", get_ip_dual_bound(optstate))
-    return optstate
+    @logmsg LogLevel(0) string("Primal bound: ", get_ip_primal_bound(outstate))
+    @logmsg LogLevel(0) string("Dual bound: ", get_ip_dual_bound(outstate))
+    return outstate, algstate
 end
 
 function optimize!(
@@ -120,7 +120,7 @@ function optimize!(
         add_lp_primal_sol!(outstate, proj_cols_on_rep(lp_primal_sol, master))
     end
 
-    return outstate
+    return outstate, algstate
 end
 
 function optimize!(
@@ -134,7 +134,7 @@ function optimize!(
     )
     algorithm = env.params.solver
     output = Algorithm.run!(algorithm, env, form, Algorithm.OptimizationInput(initstate))
-    return Algorithm.getoptstate(output)
+    return Algorithm.getoptstate(output), nothing
 end
 
 """

test/interfaces/model.jl

@@ -20,10 +20,6 @@ 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
-
 mutable struct KnapsackLibOptimizer <: BlockDecomposition.AbstractCustomOptimizer
     model::KnapsackLibModel
 end
@@ -35,6 +31,15 @@ function Coluna.Algorithm.get_units_usage(opt::KnapsackLibOptimizer, form) # for
     return units_usage
 end
 
+function _fixed_costs(model::KnapsackLibModel, form, env::Env)
+    costs = Float64[]
+    for j in 1:length(model.costs)
+        cost = Coluna.MathProg.getcurcost(form, _getvarid(model, form, env, j))
+        push!(costs, cost < 0 ? -cost : 0)
+    end
+    return costs
+end
+
 function _scale_to_int(vals...)
     return map(x -> Integer(round(10000x)), vals)
 end
@@ -45,7 +50,7 @@ function Coluna.Algorithm.run!(
     opt::KnapsackLibOptimizer, env::Coluna.Env, form::Coluna.MathProg.Formulation,
     input::Coluna.Algorithm.OptimizationInput; kw...
 )
-    costs = -[Coluna.MathProg.getcurcost(form, _getvarid(opt.model, form, env, j)) for j in 1:length(opt.model.costs)]
+    costs = _fixed_costs(opt.model, form, env)
     ws = _scale_to_int(opt.model.capacity, opt.model.weights...)
     cs = _scale_to_int(costs...)
     items = [KnapItem(w,c) for (w,c) in zip(ws[2:end], cs)]
@@ -60,8 +65,10 @@ function Coluna.Algorithm.run!(
     varvals = Float64[]
 
     for j in selected
-        push!(varids, _getvarid(opt.model, form, env, j))
-        push!(varvals, 1)
+        if costs[j] > 0
+            push!(varids, _getvarid(opt.model, form, env, j))
+            push!(varvals, 1)
+        end
     end
 
     push!(varids, setup_var_id)

test/runtests.jl

@@ -35,6 +35,7 @@ include("bound_callback_tests.jl")
 include("optimizer_with_attributes_test.jl")
 include("subproblem_solvers_tests.jl")
 include("custom_var_cuts_tests.jl")
+include("sol_disaggregation_tests.jl")
 
 rng = MersenneTwister(1234123)
 
@@ -80,3 +81,7 @@ end
 @testset "Custom Variables and Cuts" begin
     custom_var_cuts_test()
 end
+
+@testset "Solution Disaggregation" begin
+    sol_disaggregation_tests()
+end

test/sol_disaggregation_tests.jl

@@ -0,0 +1,39 @@
+function sol_disaggregation_tests()
+    I = 1:3
+    @axis(BinsType, [1])
+
+    w = [2, 5, 7]
+    Q = 8
+
+    coluna = JuMP.optimizer_with_attributes(
+        Coluna.Optimizer,
+        "params" => Coluna.Params(solver = ColumnGeneration()),
+        "default_optimizer" => GLPK.Optimizer
+    )
+
+    model = BlockModel(coluna)
+
+    @variable(model, x[k in BinsType, i in I], Bin)
+    @variable(model, y[k in BinsType], Bin)
+
+    @constraint(model, sp[i in I], sum(x[k, i] for k in BinsType) == 1)
+    @constraint(model, ks[k in BinsType], sum(w[i] * x[k, i] for i in I) - y[k] * Q <= 0)
+
+    @objective(model, Min, sum(y[k] for k in BinsType))
+
+    @dantzig_wolfe_decomposition(model, dec, BinsType)
+    subproblems = BlockDecomposition.getsubproblems(dec)
+    specify!.(subproblems, lower_multiplicity = 0, upper_multiplicity = BD.length(I)) # we use at most 3 bins
+
+    JuMP.optimize!(model)
+
+    for k in BinsType
+        bins = BD.getsolutions(model, k)
+        for bin in bins
+            @test BD.value(bin) == 1.0 # value of the master column variable
+            @test BD.value(bin, x[k, 1]) == BD.value(bin, x[k, 2]) # x[1,1] and x[1,2] in the same bin
+            @test BD.value(bin, x[k, 1]) != BD.value(bin, x[k, 3]) # only x[1,3] in its bin
+            @test BD.value(bin, x[k, 2]) != BD.value(bin, x[k, 3]) # only x[1,3] in its bin
+        end
+    end
+end