diff --git a/Project.toml b/Project.toml index 7d69f7f..ff9d8c4 100644 --- a/Project.toml +++ b/Project.toml @@ -23,10 +23,12 @@ version = "0.1.2" MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb" PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [compat] MathOptInterface = "1.34" PrecompileTools = "1" +SparseArrays = "1" Test = "1.6" julia = "1.6" diff --git a/src/MOI_wrapper.jl b/src/MOI_wrapper.jl index bf37839..d84bb24 100644 --- a/src/MOI_wrapper.jl +++ b/src/MOI_wrapper.jl @@ -18,6 +18,8 @@ # The HiGHS wrapper is released under an MIT license, a copy of which can be # found in `/thirdparty/THIRD_PARTY_LICENSES` or at https://opensource.org/licenses/MIT. +import SparseArrays + import MathOptInterface as MOI const CleverDicts = MOI.Utilities.CleverDicts @@ -621,8 +623,13 @@ end function MOI.supports( ::Optimizer, - ::Union{MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}}, -) + ::MOI.ObjectiveFunction{F}, +) where { + F<:Union{ + MOI.ScalarAffineFunction{Float64}, + MOI.ScalarQuadraticFunction{Float64}, + }, +} return true end @@ -651,6 +658,7 @@ function _check_input_data(dest::Optimizer, src::MOI.ModelLike) if attr in ( MOI.Name(), MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), + MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(), MOI.ObjectiveSense(), ) continue @@ -903,17 +911,95 @@ function _get_objective_data( objective_sense = sense == MOI.MIN_SENSE ? CUOPT_MINIMIZE : CUOPT_MAXIMIZE - objective_coefficients = zeros(Float64, numcol) F = MOI.get(src, MOI.ObjectiveFunctionType()) f_obj = MOI.get(src, MOI.ObjectiveFunction{F}()) + objective_offset, + objective_coefficients_linear, + qobj_row_offsets, + qobj_col_indices, + qobj_matrix_values = _get_objective_data(f_obj, mapping, numcol) + + return objective_sense, + objective_offset, + objective_coefficients_linear, + qobj_row_offsets, + qobj_col_indices, + qobj_matrix_values +end + +function _get_objective_data( + f_obj::MOI.ScalarAffineFunction, + mapping, + numcol::Int32, +) + objective_offset = f_obj.constant + + objective_coefficients_linear = zeros(Float64, numcol) for term in f_obj.terms - objective_coefficients[mapping[term.variable].value] += term.coefficient + i = mapping[term.variable].value + objective_coefficients_linear[i] += term.coefficient end + return objective_offset, + objective_coefficients_linear, + Int32[0], + Int32[], + Float64[] +end + +function _get_objective_data( + f_obj::MOI.ScalarQuadraticFunction, + mapping, + numcol::Int32, +) objective_offset = f_obj.constant - return objective_sense, objective_offset, objective_coefficients + objective_coefficients_linear = zeros(Float64, numcol) + for term in f_obj.affine_terms + i = mapping[term.variable].value + objective_coefficients_linear[i] += term.coefficient + end + + # Extract quadratic objective + Qtrows = Int32[] + Qtcols = Int32[] + Qtvals = Float64[] + sizehint!(Qtrows, length(f_obj.quadratic_terms)) + sizehint!(Qtcols, length(f_obj.quadratic_terms)) + sizehint!(Qtvals, length(f_obj.quadratic_terms)) + for qterm in f_obj.quadratic_terms + i = mapping[qterm.variable_1].value + j = mapping[qterm.variable_2].value + v = qterm.coefficient + # MOI stores quadratic functions as `¹/₂ xᵀQx + aᵀx + b`, with `Q` symmetric... + # (https://jump.dev/MathOptInterface.jl/stable/reference/standard_form/#MathOptInterface.ScalarQuadraticFunction) + # ... whereas cuOpt expects a QP objective of the form `¹xᵀQx + aᵀx + b`, + # where `Q` need not be symmetric + # --> we need to scale diagonal coeffs. by ¹/₂ to match cuOpt convention + if i == j + v /= 2 + end + + # We are building a COO of Qᵀ --> swap i and j + push!(Qtrows, j) + push!(Qtcols, i) + push!(Qtvals, v) + end + # CSC of Qᵀ is CSR of Q + Qt = SparseArrays.sparse(Qtrows, Qtcols, Qtvals, numcol, numcol) + qobj_matrix_values = Qt.nzval + qobj_row_offsets = Qt.colptr + qobj_col_indices = Qt.rowval + # Revert to 0-based indexing + qobj_row_offsets .-= Int32(1) + qobj_col_indices .-= Int32(1) + + return objective_offset, + objective_coefficients_linear, + qobj_row_offsets, + qobj_col_indices, + qobj_matrix_values end function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike) @@ -970,27 +1056,63 @@ function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike) has_integrality = true end - objective_sense, objective_offset, objective_coefficients = - _get_objective_data(dest, src, mapping, numcol) - - ref_problem = Ref{cuOptOptimizationProblem}() - ret = cuOptCreateRangedProblem( - numrow, - numcol, - objective_sense, - objective_offset, - objective_coefficients, - constraint_matrix_row_offsets, - constraint_matrix_column_indices, - constraint_matrix_coefficients, - rowlower, - rowupper, - collower, - colupper, - var_type, - ref_problem, - ) - _check_ret(ret, "cuOptCreateRangedProblem") + objective_sense, + objective_offset, + objective_coefficients, + qobj_row_offsets, + qobj_col_indices, + qobj_matrix_values = _get_objective_data(dest, src, mapping, numcol) + + # Is this a QP or an LP? + has_quadratic_objective = length(qobj_matrix_values) > 0 + if has_quadratic_objective && has_integrality + error( + "cuOpt does not support models with quadratic objectives _and_ integer variables", + ) + end + + if has_quadratic_objective + ref_problem = Ref{cuOptOptimizationProblem}() + ret = cuOptCreateQuadraticRangedProblem( + numrow, + numcol, + objective_sense, + objective_offset, + objective_coefficients, + qobj_row_offsets, + qobj_col_indices, + qobj_matrix_values, + constraint_matrix_row_offsets, + constraint_matrix_column_indices, + constraint_matrix_coefficients, + rowlower, + rowupper, + collower, + colupper, + ref_problem, + ) + _check_ret(ret, "cuOptCreateQuadraticRangedProblem") + else + ref_problem = Ref{cuOptOptimizationProblem}() + ret = cuOptCreateRangedProblem( + numrow, + numcol, + objective_sense, + objective_offset, + objective_coefficients, + constraint_matrix_row_offsets, + constraint_matrix_column_indices, + constraint_matrix_coefficients, + rowlower, + rowupper, + collower, + colupper, + var_type, + ref_problem, + ) + _check_ret(ret, "cuOptCreateRangedProblem") + end + dest.cuopt_problem = ref_problem[] ref_settings = Ref{cuOptSolverSettings}() diff --git a/src/cuOpt.jl b/src/cuOpt.jl index 6a3d96b..bb260d2 100644 --- a/src/cuOpt.jl +++ b/src/cuOpt.jl @@ -46,7 +46,7 @@ function __init__() error("Failed to get cuOpt library version (status code: $status)") end version = VersionNumber(major[], minor[], patch[]) - min, max = v"25.08", v"25.13" + min, max = v"25.12", v"25.13" if !(min <= version < max) error( "Incompatible cuOpt library version. Got $version, but supported versions are [$min, $max)", diff --git a/test/MOI_wrapper.jl b/test/MOI_wrapper.jl index df48a73..4136402 100644 --- a/test/MOI_wrapper.jl +++ b/test/MOI_wrapper.jl @@ -60,6 +60,10 @@ function test_runtests_cache_optimizer() "test_constraint_ZeroOne_bounds_3", # Upstream bug: https://github.com/NVIDIA/cuopt/issues/112 "test_solve_TerminationStatus_DUAL_INFEASIBLE", + # Upstream bug: https://github.com/NVIDIA/cuopt/issues/759 + # (cuOpt crashes when given a QP with no linear constraints) + "test_objective_qp_ObjectiveFunction_zero_ofdiag", + "test_objective_qp_ObjectiveFunction_edge_cases", ], ) return @@ -78,6 +82,134 @@ function test_air05() return end +function test_qp_objective() + model = MOI.instantiate(cuOpt.Optimizer; with_cache_type = Float64) + MOI.set(model, MOI.Silent(), true) + x, _ = MOI.add_constrained_variable( + model, + (MOI.GreaterThan(0.0), MOI.LessThan(1.0)), + ) + y, _ = MOI.add_constrained_variable( + model, + (MOI.GreaterThan(0.0), MOI.LessThan(1.0)), + ) + + # x + y == 1 + MOI.add_constraint( + model, + MOI.ScalarAffineFunction( + [MOI.ScalarAffineTerm(1.0, x), MOI.ScalarAffineTerm(1.0, y)], + 0.0, + ), + MOI.EqualTo(1.0), + ) + + F = MOI.ScalarQuadraticFunction{Float64} + MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE) + + # 1. Homogeneous QP with diagonal objective + # Min ¹/₂ * (x² + y²) s.t. x+y == 1 + fobj = MOI.ScalarQuadraticFunction( + [ + MOI.ScalarQuadraticTerm(1.0, x, x), + MOI.ScalarQuadraticTerm(1.0, y, y), + ], + MOI.ScalarAffineTerm{Float64}[], + 0.0, + ) + MOI.set(model, MOI.ObjectiveFunction{F}(), fobj) + MOI.optimize!(model) + + @test isapprox( + MOI.get(model, MOI.VariablePrimal(), x), + 0.5; + atol = 1e-4, + rtol = 1e-4, + ) + @test isapprox( + MOI.get(model, MOI.VariablePrimal(), y), + 0.5; + atol = 1e-4, + rtol = 1e-4, + ) + # ¹/₂ * (2 * 0.5² + 2*0.5²) == 0.5 + @test isapprox( + MOI.get(model, MOI.ObjectiveValue()), + 0.25; + atol = 1e-4, + rtol = 1e-4, + ) + + # Change diagonal coefficients + # Min ¹/₂ * (2x² + 2y²) s.t. x+y == 1 + fobj = MOI.ScalarQuadraticFunction( + [ + MOI.ScalarQuadraticTerm(2.0, x, x), + MOI.ScalarQuadraticTerm(2.0, y, y), + ], + MOI.ScalarAffineTerm{Float64}[], + 0.0, + ) + MOI.set(model, MOI.ObjectiveFunction{F}(), fobj) + MOI.optimize!(model) + # Same solution, but different objective value + # ¹/₂ * (2 * 0.5² + 2*0.5²) == 0.5 + @test isapprox( + MOI.get(model, MOI.ObjectiveValue()), + 0.5; + atol = 1e-4, + rtol = 1e-4, + ) + @test isapprox( + MOI.get(model, MOI.VariablePrimal(), x), + 0.5; + atol = 1e-4, + rtol = 1e-4, + ) + @test isapprox( + MOI.get(model, MOI.VariablePrimal(), y), + 0.5; + atol = 1e-4, + rtol = 1e-4, + ) + + # QP with off-diagonal term + # (x - y - 2)² = x² - 2xy + y² - 4x + 4y + 4 + # = ¹/₂ (2x² - xy - yx + 2y²) + (-4x + 4y) + 4 + # Solution is (x, y) = (1, 0) + fobj = MOI.ScalarQuadraticFunction( + [ + MOI.ScalarQuadraticTerm(2.0, x, x), + MOI.ScalarQuadraticTerm(2.0, y, y), + MOI.ScalarQuadraticTerm(-1.0, x, y), + ], + [MOI.ScalarAffineTerm(-4.0, x), MOI.ScalarAffineTerm(+4.0, y)], + 4.0, + ) + MOI.set(model, MOI.ObjectiveFunction{F}(), fobj) + MOI.optimize!(model) + @test isapprox( + MOI.get(model, MOI.ObjectiveValue()), + 1; + atol = 1e-4, + rtol = 1e-4, + ) + @test isapprox( + MOI.get(model, MOI.VariablePrimal(), x), + 1.0; + atol = 1e-4, + rtol = 1e-4, + ) + @test isapprox( + MOI.get(model, MOI.VariablePrimal(), y), + 0.0; + atol = 1e-4, + rtol = 1e-4, + ) + + return nothing +end + end # TestMOIWrapper TestMOIWrapper.runtests()