From 873510230551b7ff3279d547ec0d36be1c6411c0 Mon Sep 17 00:00:00 2001 From: odow Date: Wed, 6 Mar 2024 12:10:29 +1300 Subject: [PATCH 1/8] Check feasibility of solutions in PrimalStatus and DualStatus --- src/MOI_wrapper/MOI_wrapper.jl | 34 ++++++++++++++++++++++++++-------- 1 file changed, 26 insertions(+), 8 deletions(-) diff --git a/src/MOI_wrapper/MOI_wrapper.jl b/src/MOI_wrapper/MOI_wrapper.jl index 313d8e0..092fecf 100644 --- a/src/MOI_wrapper/MOI_wrapper.jl +++ b/src/MOI_wrapper/MOI_wrapper.jl @@ -2831,8 +2831,19 @@ function MOI.get(model::Optimizer, attr::MOI.PrimalStatus) valueP = Ref{Cint}() ret = GRBgetintattr(model, "SolCount", valueP) _check_ret(model, ret) - return 1 <= attr.result_index <= valueP[] ? MOI.FEASIBLE_POINT : - MOI.NO_SOLUTION + if !(1 <= attr.result_index <= valueP[]) + return MOI.NO_SOLUTION + end + doubleP = Ref{Cdouble}() + ret = GRBgetdblattr(model, "MaxVio", doubleP) + _check_ret(model, ret) + if doubleP[] < 1e-8 + return MOI.FEASIBLE_POINT + elseif doubleP[] < 1e-6 + return MOI.NEARLY_FEASIBLE_POINT + else + return MOI.INFEASIBLE_POINT + end end function _has_dual_ray(model::Optimizer) @@ -2869,18 +2880,25 @@ function MOI.get(model::Optimizer, attr::MOI.DualStatus) valueP = Ref{Cint}() ret = GRBgetintattr(model, "SolCount", valueP) _check_ret(model, ret) - if valueP[] == 0 + if !(1 <= attr.result_index <= valueP[]) return MOI.NO_SOLUTION end # But unfortunately, even if SolCount is 1, sometimes a dual solution is not # available. The only way to check this is to query directly. - valueP = Ref{Cdouble}() - ret = GRBgetdblattrelement(model, "RC", 0, valueP) - if ret == 0 - return MOI.FEASIBLE_POINT - else # Something went wrong + doubleP = Ref{Cdouble}() + ret = GRBgetdblattrelement(model, "RC", 0, doubleP) + if ret != 0 # Something went wrong return MOI.NO_SOLUTION end + ret = GRBgetdblattr(model, "MaxVio", doubleP) + _check_ret(model, ret) + if doubleP[] < 1e-8 + return MOI.FEASIBLE_POINT + elseif doubleP[] < 1e-6 + return MOI.NEARLY_FEASIBLE_POINT + else + return MOI.INFEASIBLE_POINT + end end function _has_primal_ray(model::Optimizer) From d6b08a0805d78645ceaf0993f0e6011ecefdaa8a Mon Sep 17 00:00:00 2001 From: odow Date: Wed, 6 Mar 2024 13:14:27 +1300 Subject: [PATCH 2/8] Update --- src/MOI_wrapper/MOI_wrapper.jl | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/MOI_wrapper/MOI_wrapper.jl b/src/MOI_wrapper/MOI_wrapper.jl index 092fecf..f2a0e22 100644 --- a/src/MOI_wrapper/MOI_wrapper.jl +++ b/src/MOI_wrapper/MOI_wrapper.jl @@ -2834,6 +2834,10 @@ function MOI.get(model::Optimizer, attr::MOI.PrimalStatus) if !(1 <= attr.result_index <= valueP[]) return MOI.NO_SOLUTION end + if term in (MOI.OPTIMAL, MOI.SOLUTION_LIMIT) + return MOI.FEASIBLE_POINT + end + # Feasibility of solution is unknown. Check for violations. doubleP = Ref{Cdouble}() ret = GRBgetdblattr(model, "MaxVio", doubleP) _check_ret(model, ret) @@ -2890,6 +2894,10 @@ function MOI.get(model::Optimizer, attr::MOI.DualStatus) if ret != 0 # Something went wrong return MOI.NO_SOLUTION end + if term in (MOI.OPTIMAL, MOI.SOLUTION_LIMIT) + return MOI.FEASIBLE_POINT + end + # Feasibility of solution is unknown. Check for violations. ret = GRBgetdblattr(model, "MaxVio", doubleP) _check_ret(model, ret) if doubleP[] < 1e-8 From ad7620de4249f89c0664e2ecd80373945371b028 Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Wed, 6 Mar 2024 16:21:40 +1300 Subject: [PATCH 3/8] Update MOI_wrapper.jl --- src/MOI_wrapper/MOI_wrapper.jl | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) diff --git a/src/MOI_wrapper/MOI_wrapper.jl b/src/MOI_wrapper/MOI_wrapper.jl index f2a0e22..9f02941 100644 --- a/src/MOI_wrapper/MOI_wrapper.jl +++ b/src/MOI_wrapper/MOI_wrapper.jl @@ -2843,11 +2843,8 @@ function MOI.get(model::Optimizer, attr::MOI.PrimalStatus) _check_ret(model, ret) if doubleP[] < 1e-8 return MOI.FEASIBLE_POINT - elseif doubleP[] < 1e-6 - return MOI.NEARLY_FEASIBLE_POINT - else - return MOI.INFEASIBLE_POINT end + return MOI.UNKNOWN_RESULT_STATUS end function _has_dual_ray(model::Optimizer) @@ -2902,11 +2899,8 @@ function MOI.get(model::Optimizer, attr::MOI.DualStatus) _check_ret(model, ret) if doubleP[] < 1e-8 return MOI.FEASIBLE_POINT - elseif doubleP[] < 1e-6 - return MOI.NEARLY_FEASIBLE_POINT - else - return MOI.INFEASIBLE_POINT end + return MOI.UNKNOWN_RESULT_STATUS end function _has_primal_ray(model::Optimizer) From 286341e384d4fe250f6d13ed326d5b11dbbde1fb Mon Sep 17 00:00:00 2001 From: odow Date: Fri, 8 Mar 2024 11:35:59 +1300 Subject: [PATCH 4/8] Add _is_primal_feasible_to_tolerance --- src/MOI_wrapper/MOI_wrapper.jl | 64 ++++++++++++++++++++++++++-------- 1 file changed, 50 insertions(+), 14 deletions(-) diff --git a/src/MOI_wrapper/MOI_wrapper.jl b/src/MOI_wrapper/MOI_wrapper.jl index 9f02941..0daced5 100644 --- a/src/MOI_wrapper/MOI_wrapper.jl +++ b/src/MOI_wrapper/MOI_wrapper.jl @@ -2833,20 +2833,60 @@ function MOI.get(model::Optimizer, attr::MOI.PrimalStatus) _check_ret(model, ret) if !(1 <= attr.result_index <= valueP[]) return MOI.NO_SOLUTION - end - if term in (MOI.OPTIMAL, MOI.SOLUTION_LIMIT) + elseif term in (MOI.OPTIMAL, MOI.SOLUTION_LIMIT) + # Assume that if Gurobi tells us the problem is optimal, or if it found + # too manny solutions, then we have a feasible point. return MOI.FEASIBLE_POINT - end - # Feasibility of solution is unknown. Check for violations. - doubleP = Ref{Cdouble}() - ret = GRBgetdblattr(model, "MaxVio", doubleP) - _check_ret(model, ret) - if doubleP[] < 1e-8 + elseif _is_primal_feasible_to_tolerance(model) + # Feasibility of solution is unknown. Check for violations. return MOI.FEASIBLE_POINT end return MOI.UNKNOWN_RESULT_STATUS end +""" + _is_primal_feasible_to_tolerance(model::Optimizer) + +This function checks whether a solution stored in `model` is feasible according +to: + + 1. `max(BoundVio, ConstrVio) <= FeasibilityTol` + 2. `IntVio <= IntFeasTol` (if applicable) + +For a discussion on why this is needed, see: + + * https://github.com/jump-dev/Gurobi.jl/pull/546 + * https://github.com/jump-dev/Gurobi.jl/pull/548 + +Ideally, a future version of Gurobi would provide a C function for this. +""" +function _is_primal_feasible_to_tolerance(model::Optimizer) + boundVioP, constrVioP = Ref{Cdouble}(0.0), Ref{Cdouble}(0.0) + ret = GRBgetdlbattr(model, "ConstrVio", constrVioP) + _check_ret(model, ret) + ret = GRBgetdlbattr(model, "BoundVio", boundVioP) + _check_ret(model, ret) + feasibilityTolP = Ref{Cdouble}(0.0) + ret = GRBgetdlbparam(model, "FeasibilityTol", feasibilityTolP) + _check_ret(model, ret) + if max(boundVioP[], constrVioP[]) > feasibilityTolP[] + return false + end + intVioP = Ref{Cdouble}(0.0) + ret = GRBgetdlbattr(model, "IntVio", intVioP) + if ret == GRB_ERROR_DATA_NOT_AVAILABLE + return true # IntVio is available only for Int models + end + _check_ret(model, ret) + intFeasTolP = Ref{Cdouble}(0.0) + ret = GRBgetdlbparam(model, "IntFeasTol", intFeasTolP) + _check_ret(model, ret) + if intVioP[] > intFeasTolP[] + return false + end + return true +end + function _has_dual_ray(model::Optimizer) # Note: for performance reasons, we try to get 1 element because for # some versions of Gurobi, we cannot query 0 elements without error. @@ -2894,12 +2934,8 @@ function MOI.get(model::Optimizer, attr::MOI.DualStatus) if term in (MOI.OPTIMAL, MOI.SOLUTION_LIMIT) return MOI.FEASIBLE_POINT end - # Feasibility of solution is unknown. Check for violations. - ret = GRBgetdblattr(model, "MaxVio", doubleP) - _check_ret(model, ret) - if doubleP[] < 1e-8 - return MOI.FEASIBLE_POINT - end + # There is no easy equivalence to _is_primal_feasible_to_tolerance, so we + # punt and return UNKNOWN_RESULT_STATUS. return MOI.UNKNOWN_RESULT_STATUS end From 8893511f8b9f1f6d3579a6f3c1b3cf98ff64b8e5 Mon Sep 17 00:00:00 2001 From: odow Date: Fri, 8 Mar 2024 11:38:09 +1300 Subject: [PATCH 5/8] Fix spelling --- src/MOI_wrapper/MOI_wrapper.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/MOI_wrapper/MOI_wrapper.jl b/src/MOI_wrapper/MOI_wrapper.jl index 0daced5..4e0fb0e 100644 --- a/src/MOI_wrapper/MOI_wrapper.jl +++ b/src/MOI_wrapper/MOI_wrapper.jl @@ -2862,24 +2862,24 @@ Ideally, a future version of Gurobi would provide a C function for this. """ function _is_primal_feasible_to_tolerance(model::Optimizer) boundVioP, constrVioP = Ref{Cdouble}(0.0), Ref{Cdouble}(0.0) - ret = GRBgetdlbattr(model, "ConstrVio", constrVioP) + ret = GRBgetdblattr(model, "ConstrVio", constrVioP) _check_ret(model, ret) - ret = GRBgetdlbattr(model, "BoundVio", boundVioP) + ret = GRBgetdblattr(model, "BoundVio", boundVioP) _check_ret(model, ret) feasibilityTolP = Ref{Cdouble}(0.0) - ret = GRBgetdlbparam(model, "FeasibilityTol", feasibilityTolP) + ret = GRBgetdblparam(model, "FeasibilityTol", feasibilityTolP) _check_ret(model, ret) if max(boundVioP[], constrVioP[]) > feasibilityTolP[] return false end intVioP = Ref{Cdouble}(0.0) - ret = GRBgetdlbattr(model, "IntVio", intVioP) + ret = GRBgetdblattr(model, "IntVio", intVioP) if ret == GRB_ERROR_DATA_NOT_AVAILABLE return true # IntVio is available only for Int models end _check_ret(model, ret) intFeasTolP = Ref{Cdouble}(0.0) - ret = GRBgetdlbparam(model, "IntFeasTol", intFeasTolP) + ret = GRBgetdblparam(model, "IntFeasTol", intFeasTolP) _check_ret(model, ret) if intVioP[] > intFeasTolP[] return false From 4bfe9f14b56cfdbae3a35e29a049558ee8cdec2d Mon Sep 17 00:00:00 2001 From: odow Date: Fri, 8 Mar 2024 12:14:57 +1300 Subject: [PATCH 6/8] Add tests --- src/MOI_wrapper/MOI_wrapper.jl | 23 ++++++++------- test/MOI/MOI_wrapper.jl | 51 ++++++++++++++++++++++++++++++++++ 2 files changed, 62 insertions(+), 12 deletions(-) diff --git a/src/MOI_wrapper/MOI_wrapper.jl b/src/MOI_wrapper/MOI_wrapper.jl index 4e0fb0e..9ea0f95 100644 --- a/src/MOI_wrapper/MOI_wrapper.jl +++ b/src/MOI_wrapper/MOI_wrapper.jl @@ -2863,28 +2863,27 @@ Ideally, a future version of Gurobi would provide a C function for this. function _is_primal_feasible_to_tolerance(model::Optimizer) boundVioP, constrVioP = Ref{Cdouble}(0.0), Ref{Cdouble}(0.0) ret = GRBgetdblattr(model, "ConstrVio", constrVioP) + if ret == GRB_ERROR_DATA_NOT_AVAILABLE + return false # No solution present + end _check_ret(model, ret) ret = GRBgetdblattr(model, "BoundVio", boundVioP) _check_ret(model, ret) - feasibilityTolP = Ref{Cdouble}(0.0) - ret = GRBgetdblparam(model, "FeasibilityTol", feasibilityTolP) - _check_ret(model, ret) - if max(boundVioP[], constrVioP[]) > feasibilityTolP[] - return false - end intVioP = Ref{Cdouble}(0.0) ret = GRBgetdblattr(model, "IntVio", intVioP) if ret == GRB_ERROR_DATA_NOT_AVAILABLE - return true # IntVio is available only for Int models + intVioP[] = 0.0 + else + _check_ret(model, ret) end + feasibilityTolP = Ref{Cdouble}(0.0) + ret = GRBgetdblparam(model.env, "FeasibilityTol", feasibilityTolP) _check_ret(model, ret) intFeasTolP = Ref{Cdouble}(0.0) - ret = GRBgetdblparam(model, "IntFeasTol", intFeasTolP) + ret = GRBgetdblparam(model.env, "IntFeasTol", intFeasTolP) _check_ret(model, ret) - if intVioP[] > intFeasTolP[] - return false - end - return true + return max(boundVioP[], constrVioP[]) <= feasibilityTolP[] && + intVioP[] <= intFeasTolP[] end function _has_dual_ray(model::Optimizer) diff --git a/test/MOI/MOI_wrapper.jl b/test/MOI/MOI_wrapper.jl index 4d01288..5856b64 100644 --- a/test/MOI/MOI_wrapper.jl +++ b/test/MOI/MOI_wrapper.jl @@ -898,6 +898,57 @@ function test_delete_indicator() return end +function test_is_primal_feasible_to_tolerance() + model = Gurobi.Optimizer(GRB_ENV) + MOI.Utilities.loadfromstring!( + model, + """ + variables: x, z + maxobjective: 2.0x + 1000.0 + c1: x + -1_333_333.12345 * z <= 0.0 + c2: z <= 0.000001 + c3: z in ZeroOne() + """, + ) + MOI.optimize!(model) + @test Gurobi._is_primal_feasible_to_tolerance(model) + return +end + +function test_is_primal_feasible_to_tolerance_infeasible() + model = Gurobi.Optimizer(GRB_ENV) + MOI.Utilities.loadfromstring!( + model, + """ + variables: x, z + maxobjective: 2.0x + 1000.0 + x + -1.0 * z <= 0.0 + z <= 0.5 + x >= 0.5 + z in ZeroOne() + """, + ) + MOI.optimize!(model) + @test !Gurobi._is_primal_feasible_to_tolerance(model) + return end +function test_is_primal_feasible_to_tolerance_lp() + model = Gurobi.Optimizer(GRB_ENV) + MOI.Utilities.loadfromstring!( + model, + """ + variables: x + minobjective: x + c1: x >= 0.0 + c2: 2x >= 1.0 + """, + ) + MOI.optimize!(model) + @test Gurobi._is_primal_feasible_to_tolerance(model) + return +end + +end # TestMOIWrapper + TestMOIWrapper.runtests() From cd4c2701284b6d70be9798c6170c52c45295afcf Mon Sep 17 00:00:00 2001 From: odow Date: Fri, 8 Mar 2024 15:27:00 +1300 Subject: [PATCH 7/8] Add more tests --- test/MOI/MOI_wrapper.jl | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) diff --git a/test/MOI/MOI_wrapper.jl b/test/MOI/MOI_wrapper.jl index 5856b64..e8d2742 100644 --- a/test/MOI/MOI_wrapper.jl +++ b/test/MOI/MOI_wrapper.jl @@ -949,6 +949,34 @@ function test_is_primal_feasible_to_tolerance_lp() return end +function test_primal_feasible_status() + model = Gurobi.Optimizer(GRB_ENV) + MOI.set(model, MOI.Silent(), true) + MOI.set(model, MOI.RawOptimizerAttribute("Heuristics"), 0) + MOI.set(model, MOI.RawOptimizerAttribute("Presolve"), 0) + MOI.set(model, MOI.RawOptimizerAttribute("IterationLimit"), 0) + N = 100 + x = MOI.add_variables(model, N) + MOI.add_constraint.(model, x, MOI.ZeroOne()) + MOI.set.(model, MOI.VariablePrimalStart(), x, 0.0) + Random.seed!(1) + item_weights, item_values = rand(N), rand(N) + MOI.add_constraint( + model, + MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(item_weights, x), 0.0), + MOI.LessThan(10.0), + ) + MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE) + MOI.set( + model, + MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), + MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(item_values, x), 0.0), + ) + MOI.optimize!(model) + @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT + return +end + end # TestMOIWrapper TestMOIWrapper.runtests() From ae1e4a99a38b48f9bf769b73260441347640c0e2 Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Mon, 11 Mar 2024 20:13:43 +1300 Subject: [PATCH 8/8] Apply suggestions from code review --- src/MOI_wrapper/MOI_wrapper.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/MOI_wrapper/MOI_wrapper.jl b/src/MOI_wrapper/MOI_wrapper.jl index 9ea0f95..94a49b1 100644 --- a/src/MOI_wrapper/MOI_wrapper.jl +++ b/src/MOI_wrapper/MOI_wrapper.jl @@ -2925,8 +2925,8 @@ function MOI.get(model::Optimizer, attr::MOI.DualStatus) end # But unfortunately, even if SolCount is 1, sometimes a dual solution is not # available. The only way to check this is to query directly. - doubleP = Ref{Cdouble}() - ret = GRBgetdblattrelement(model, "RC", 0, doubleP) + valueP = Ref{Cdouble}() + ret = GRBgetdblattrelement(model, "RC", 0, valueP) if ret != 0 # Something went wrong return MOI.NO_SOLUTION end