Enable algorithms to update the ideal_point during a solve (#102)

Author: odow

src/MultiObjectiveAlgorithms.jl

@@ -560,15 +560,13 @@ function MOI.delete(model::Optimizer, ci::MOI.ConstraintIndex)
 end
 
 function _compute_ideal_point(model::Optimizer, start_time)
-    objectives = MOI.Utilities.eachscalar(model.f)
-    model.ideal_point = fill(NaN, length(objectives))
-    if !MOI.get(model, ComputeIdealPoint())
-        return
-    end
-    for (i, f) in enumerate(objectives)
+    for (i, f) in enumerate(MOI.Utilities.eachscalar(model.f))
         if _time_limit_exceeded(model, start_time)
             return
         end
+        if !isnan(model.ideal_point[i])
+            continue  # The algorithm already updated this information
+        end
         MOI.set(model.inner, MOI.ObjectiveFunction{typeof(f)}(), f)
         MOI.optimize!(model.inner)
         status = MOI.get(model.inner, MOI.TerminationStatus())
@@ -585,15 +583,21 @@ function MOI.optimize!(model::Optimizer)
     model.termination_status = MOI.OPTIMIZE_NOT_CALLED
     if model.f === nothing
         model.termination_status = MOI.INVALID_MODEL
+        empty!(model.ideal_point)
         return
     end
+    # We need to clear the ideal point prior to starting the solve. Algorithms
+    # may update this during the solve, otherwise we will update it at the end.
+    model.ideal_point = fill(NaN, MOI.output_dimension(model.f))
     algorithm = something(model.algorithm, default(Algorithm()))
     status, solutions = optimize_multiobjective!(algorithm, model)
     model.termination_status = status
-    _compute_ideal_point(model, start_time)
     if solutions !== nothing
         model.solutions = solutions
     end
+    if MOI.get(model, ComputeIdealPoint())
+        _compute_ideal_point(model, start_time)
+    end
     if MOI.supports(model.inner, MOI.TimeLimitSec())
         MOI.set(model.inner, MOI.TimeLimitSec(), nothing)
     end

src/algorithms/Dichotomy.jl

@@ -54,10 +54,6 @@ function MOI.get(alg::Dichotomy, attr::SolutionLimit)
     return something(alg.solution_limit, default(alg, attr))
 end
 
-function _solve_weighted_sum(model::Optimizer, alg::Dichotomy, weight::Float64)
-    return _solve_weighted_sum(model, alg, [weight, 1 - weight])
-end
-
 function _solve_weighted_sum(
     model::Optimizer,
     ::Dichotomy,
@@ -88,15 +84,17 @@ function optimize_multiobjective!(algorithm::Dichotomy, model::Optimizer)
         return status, [solution]
     end
     solutions = Dict{Float64,SolutionPoint}()
-    for w in (0.0, 1.0)
+    for (i, w) in (1 => 1.0, 2 => 0.0)
         if _time_limit_exceeded(model, start_time)
             return MOI.TIME_LIMIT, nothing
         end
-        status, solution = _solve_weighted_sum(model, algorithm, w)
+        status, solution = _solve_weighted_sum(model, algorithm, [w, 1.0 - w])
         if !_is_scalar_status_optimal(status)
             return status, nothing
         end
         solutions[w] = solution
+        # We already have enough information here to update the ideal point.
+        model.ideal_point[i] = solution.y[i]
     end
     queue = Tuple{Float64,Float64}[]
     if !(solutions[0.0] ≈ solutions[1.0])
@@ -112,7 +110,7 @@ function optimize_multiobjective!(algorithm::Dichotomy, model::Optimizer)
         (a, b) = popfirst!(queue)
         y_d = solutions[a].y .- solutions[b].y
         w = y_d[2] / (y_d[2] - y_d[1])
-        status, solution = _solve_weighted_sum(model, algorithm, w)
+        status, solution = _solve_weighted_sum(model, algorithm, [w, 1.0 - w])
         if !_is_scalar_status_optimal(status)
             # Exit the solve with some error.
             return status, nothing

src/algorithms/DominguezRios.jl

@@ -155,6 +155,7 @@ function optimize_multiobjective!(algorithm::DominguezRios, model::Optimizer)
         if solutions !== nothing
             solutions = [SolutionPoint(s.x, -s.y) for s in solutions]
         end
+        model.ideal_point .*= -1
         return status, solutions
     end
     n = MOI.output_dimension(model.f)
@@ -173,6 +174,7 @@ function optimize_multiobjective!(algorithm::DominguezRios, model::Optimizer)
         end
         _, Y = _compute_point(model, variables, f_i)
         yI[i] = Y
+        model.ideal_point[i] = Y
         rev_sense = sense == MOI.MIN_SENSE ? MOI.MAX_SENSE : MOI.MIN_SENSE
         MOI.set(model.inner, MOI.ObjectiveSense(), rev_sense)
         MOI.optimize!(model.inner)

src/algorithms/EpsilonConstraint.jl

@@ -73,7 +73,7 @@ function optimize_multiobjective!(
     if MOI.output_dimension(model.f) != 2
         error("EpsilonConstraint requires exactly two objectives")
     end
-    # Compute the bounding box ofthe objectives using Hierarchical().
+    # Compute the bounding box of the objectives using Hierarchical().
     alg = Hierarchical()
     MOI.set.(Ref(alg), ObjectivePriority.(1:2), [1, 0])
     status, solution_1 = optimize_multiobjective!(alg, model)
@@ -87,22 +87,27 @@ function optimize_multiobjective!(
     end
     a, b = solution_1[1].y[1], solution_2[1].y[1]
     left, right = min(a, b), max(a, b)
+    sense = MOI.get(model.inner, MOI.ObjectiveSense())
+    if sense == MOI.MIN_SENSE
+        model.ideal_point .= min.(solution_1[1].y, solution_2[1].y)
+    else
+        model.ideal_point .= max.(solution_1[1].y, solution_2[1].y)
+    end
     # Compute the epsilon that we will be incrementing by each iteration
     ε = MOI.get(algorithm, EpsilonConstraintStep())
     n_points = MOI.get(algorithm, SolutionLimit())
     if n_points != default(algorithm, SolutionLimit())
         ε = abs(right - left) / (n_points - 1)
     end
-    solutions = SolutionPoint[]
+    solutions = SolutionPoint[only(solution_1), only(solution_2)]
     f1, f2 = MOI.Utilities.eachscalar(model.f)
     MOI.set(model.inner, MOI.ObjectiveFunction{typeof(f2)}(), f2)
     # Add epsilon constraint
-    sense = MOI.get(model.inner, MOI.ObjectiveSense())
     variables = MOI.get(model.inner, MOI.ListOfVariableIndices())
     SetType, bound = if sense == MOI.MIN_SENSE
-        MOI.LessThan{Float64}, right
+        MOI.LessThan{Float64}, right - ε
     else
-        MOI.GreaterThan{Float64}, left
+        MOI.GreaterThan{Float64}, left + ε
     end
     constant = MOI.constant(f1, Float64)
     ci = MOI.Utilities.normalize_and_add_constraint(
@@ -113,7 +118,7 @@ function optimize_multiobjective!(
     )
     bound -= constant
     status = MOI.OPTIMAL
-    for _ in 1:n_points
+    for _ in 3:n_points
         if _time_limit_exceeded(model, start_time)
             status = MOI.TIME_LIMIT
             break

src/algorithms/KirlikSayin.jl

@@ -87,6 +87,7 @@ function optimize_multiobjective!(algorithm::KirlikSayin, model::Optimizer)
         if solutions !== nothing
             solutions = [SolutionPoint(s.x, -s.y) for s in solutions]
         end
+        model.ideal_point .*= -1
         return status, solutions
     end
     solutions = SolutionPoint[]
@@ -111,6 +112,7 @@ function optimize_multiobjective!(algorithm::KirlikSayin, model::Optimizer)
         end
         _, Y = _compute_point(model, variables, f_i)
         yI[i] = Y + 1
+        model.ideal_point[i] = Y
         MOI.set(
             model.inner,
             MOI.ObjectiveSense(),

src/algorithms/TambyVanderpooten.jl

@@ -92,6 +92,7 @@ function optimize_multiobjective!(
         if solutions !== nothing
             solutions = [SolutionPoint(s.x, -s.y) for s in solutions]
         end
+        model.ideal_point .*= -1
         return status, solutions
     end
     warm_start_supported = false
@@ -114,6 +115,7 @@ function optimize_multiobjective!(
         end
         _, Y = _compute_point(model, variables, f_i)
         yI[i] = Y + 1
+        model.ideal_point[i] = Y
         MOI.set(model.inner, MOI.ObjectiveSense(), MOI.MAX_SENSE)
         MOI.optimize!(model.inner)
         status = MOI.get(model.inner, MOI.TerminationStatus())

test/algorithms/Chalmet.jl

@@ -68,6 +68,7 @@ function test_knapsack_min()
     @test isapprox(x_sol, X_E'; atol = 1e-6)
     y_sol = hcat([MOI.get(model, MOI.ObjectiveValue(i)) for i in 1:N]...)
     @test isapprox(y_sol, Y_N'; atol = 1e-6)
+    @test MOI.get(model, MOI.ObjectiveBound()) ≈ vec(minimum(Y_N; dims = 1))
     return
 end
 
@@ -117,6 +118,7 @@ function test_knapsack_max()
     @test isapprox(x_sol, X_E'; atol = 1e-6)
     y_sol = hcat([MOI.get(model, MOI.ObjectiveValue(i)) for i in 1:N]...)
     @test isapprox(y_sol, Y_N'; atol = 1e-6)
+    @test MOI.get(model, MOI.ObjectiveBound()) ≈ vec(maximum(Y_N; dims = 1))
     return
 end
 
@@ -195,14 +197,17 @@ function test_vector_of_variables_objective()
     end
     MOI.set(model, MOA.Algorithm(), MOA.Chalmet())
     MOI.set(model, MOI.Silent(), true)
+    MOI.set(model, MOA.ComputeIdealPoint(), false)
     x = MOI.add_variables(model, 2)
     MOI.add_constraint.(model, x, MOI.ZeroOne())
     f = MOI.VectorOfVariables(x)
     MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
     MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
     MOI.add_constraint(model, sum(1.0 * xi for xi in x), MOI.GreaterThan(1.0))
     MOI.optimize!(model)
-    MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
+    @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
+    ideal_point = MOI.get(model, MOI.ObjectiveBound())
+    @test length(ideal_point) == 2 && all(isnan, ideal_point)
     return
 end
 

test/algorithms/DominguezRios.jl

@@ -82,6 +82,7 @@ function test_knapsack_min_p3()
     @test isapprox(sort(x_sol; dims = 1), sort(X_E'; dims = 1); atol = 1e-6)
     y_sol = vcat([MOI.get(model, MOI.ObjectiveValue(i))' for i in 1:N]...)
     @test isapprox(sort(y_sol; dims = 1), sort(Y_N; dims = 1); atol = 1e-6)
+    @test MOI.get(model, MOI.ObjectiveBound()) ≈ vec(minimum(Y_N; dims = 1))
     return
 end
 
@@ -141,6 +142,7 @@ function test_knapsack_max_p3()
     @test isapprox(sort(x_sol; dims = 1), sort(X_E'; dims = 1); atol = 1e-6)
     y_sol = vcat([MOI.get(model, MOI.ObjectiveValue(i))' for i in 1:N]...)
     @test isapprox(sort(y_sol; dims = 1), sort(Y_N; dims = 1); atol = 1e-6)
+    @test MOI.get(model, MOI.ObjectiveBound()) ≈ vec(maximum(Y_N; dims = 1))
     return
 end
 
@@ -207,6 +209,7 @@ function test_knapsack_min_p4()
     @test isapprox(sort(x_sol; dims = 1), sort(X_E'; dims = 1); atol = 1e-6)
     y_sol = vcat([MOI.get(model, MOI.ObjectiveValue(i))' for i in 1:N]...)
     @test isapprox(sort(y_sol; dims = 1), sort(Y_N; dims = 1); atol = 1e-6)
+    @test MOI.get(model, MOI.ObjectiveBound()) ≈ vec(minimum(Y_N; dims = 1))
     return
 end
 
@@ -273,6 +276,7 @@ function test_knapsack_max_p4()
     @test isapprox(sort(x_sol; dims = 1), sort(X_E'; dims = 1); atol = 1e-6)
     y_sol = vcat([MOI.get(model, MOI.ObjectiveValue(i))' for i in 1:N]...)
     @test isapprox(sort(y_sol; dims = 1), sort(Y_N; dims = 1); atol = 1e-6)
+    @test MOI.get(model, MOI.ObjectiveBound()) ≈ vec(maximum(Y_N; dims = 1))
     return
 end
 

test/algorithms/EpsilonConstraint.jl

@@ -65,6 +65,7 @@ function test_biobjective_knapsack()
         Y = MOI.get(model, MOI.ObjectiveValue(i))
         @test results[round.(Int, Y)] == X
     end
+    @test MOI.get(model, MOI.ObjectiveBound()) == [956.0, 983.0]
     return
 end
 
@@ -108,6 +109,7 @@ function test_biobjective_knapsack_atol()
         Y = MOI.get(model, MOI.ObjectiveValue(i))
         @test results[round.(Int, Y)] == X
     end
+    @test MOI.get(model, MOI.ObjectiveBound()) == [955.0, 983.0]
     return
 end
 
@@ -136,11 +138,12 @@ function test_biobjective_knapsack_atol_large()
     )
     MOI.optimize!(model)
     results = Dict(
+        [955, 906] => [2, 3, 5, 6, 9, 10, 11, 14, 15, 16, 17],
         [948, 939] => [1, 2, 3, 5, 6, 8, 10, 11, 15, 16, 17],
         [934, 971] => [2, 3, 5, 6, 8, 10, 11, 12, 15, 16, 17],
         [918, 983] => [2, 3, 4, 5, 6, 8, 10, 11, 12, 16, 17],
     )
-    @test MOI.get(model, MOI.ResultCount()) == 3
+    @test MOI.get(model, MOI.ResultCount()) == 4
     for i in 1:MOI.get(model, MOI.ResultCount())
         x_sol = MOI.get(model, MOI.VariablePrimal(i), x)
         X = findall(elt -> elt > 0.9, x_sol)
@@ -191,6 +194,7 @@ function test_biobjective_knapsack_min()
         Y = MOI.get(model, MOI.ObjectiveValue(i))
         @test results[-round.(Int, Y)] == X
     end
+    @test MOI.get(model, MOI.ObjectiveBound()) == [-955.0, -983.0]
     return
 end
 
@@ -219,16 +223,18 @@ function test_biobjective_knapsack_min_solution_limit()
     )
     MOI.optimize!(model)
     results = Dict(
+        [955, 906] => [2, 3, 5, 6, 9, 10, 11, 14, 15, 16, 17],
         [943, 940] => [2, 3, 5, 6, 8, 9, 10, 11, 15, 16, 17],
         [918, 983] => [2, 3, 4, 5, 6, 8, 10, 11, 12, 16, 17],
     )
-    @test MOI.get(model, MOI.ResultCount()) == 2
+    @test MOI.get(model, MOI.ResultCount()) == 3
     for i in 1:MOI.get(model, MOI.ResultCount())
         x_sol = MOI.get(model, MOI.VariablePrimal(i), x)
         X = findall(elt -> elt > 0.9, x_sol)
         Y = MOI.get(model, MOI.ObjectiveValue(i))
         @test results[round.(Int, Y)] == X
     end
+    @test MOI.get(model, MOI.ObjectiveBound()) == [955.0, 983.0]
     return
 end
 
@@ -419,7 +425,8 @@ function test_time_limit()
     )
     MOI.optimize!(model)
     @test MOI.get(model, MOI.TerminationStatus()) == MOI.TIME_LIMIT
-    @test MOI.get(model, MOI.ResultCount()) == 0
+    # Check time limits in subsolves
+    @test_broken MOI.get(model, MOI.ResultCount()) == 0
     return
 end
 

test/algorithms/KirlikSayin.jl

@@ -79,6 +79,7 @@ function test_knapsack_min_p3()
     @test isapprox(sort(x_sol; dims = 1), sort(X_E'; dims = 1); atol = 1e-6)
     y_sol = vcat([MOI.get(model, MOI.ObjectiveValue(i))' for i in 1:N]...)
     @test isapprox(sort(y_sol; dims = 1), sort(Y_N; dims = 1); atol = 1e-6)
+    @test MOI.get(model, MOI.ObjectiveBound()) ≈ vec(minimum(Y_N; dims = 1))
     return
 end
 
@@ -138,6 +139,7 @@ function test_knapsack_max_p3()
     @test isapprox(sort(x_sol; dims = 1), sort(X_E'; dims = 1); atol = 1e-6)
     y_sol = vcat([MOI.get(model, MOI.ObjectiveValue(i))' for i in 1:N]...)
     @test isapprox(sort(y_sol; dims = 1), sort(Y_N; dims = 1); atol = 1e-6)
+    @test MOI.get(model, MOI.ObjectiveBound()) ≈ vec(maximum(Y_N; dims = 1))
     return
 end
 
@@ -204,6 +206,7 @@ function test_knapsack_min_p4()
     @test isapprox(sort(x_sol; dims = 1), sort(X_E'; dims = 1); atol = 1e-6)
     y_sol = vcat([MOI.get(model, MOI.ObjectiveValue(i))' for i in 1:N]...)
     @test isapprox(sort(y_sol; dims = 1), sort(Y_N; dims = 1); atol = 1e-6)
+    @test MOI.get(model, MOI.ObjectiveBound()) ≈ vec(minimum(Y_N; dims = 1))
     return
 end
 
@@ -270,6 +273,7 @@ function test_knapsack_max_p4()
     @test isapprox(sort(x_sol; dims = 1), sort(X_E'; dims = 1); atol = 1e-6)
     y_sol = vcat([MOI.get(model, MOI.ObjectiveValue(i))' for i in 1:N]...)
     @test isapprox(sort(y_sol; dims = 1), sort(Y_N; dims = 1); atol = 1e-6)
+    @test MOI.get(model, MOI.ObjectiveBound()) ≈ vec(maximum(Y_N; dims = 1))
     return
 end