diff --git a/conda/environments/all_cuda-129_arch-aarch64.yaml b/conda/environments/all_cuda-129_arch-aarch64.yaml index 37977ea28..6a7ffe948 100644 --- a/conda/environments/all_cuda-129_arch-aarch64.yaml +++ b/conda/environments/all_cuda-129_arch-aarch64.yaml @@ -61,6 +61,7 @@ dependencies: - requests - rmm==26.2.*,>=0.0.0a0 - scikit-build-core>=0.11.0 +- scipy>=1.13.0 - sphinx - sphinx-copybutton - sphinx-design diff --git a/conda/environments/all_cuda-129_arch-x86_64.yaml b/conda/environments/all_cuda-129_arch-x86_64.yaml index 0eaa7000a..7b4853e96 100644 --- a/conda/environments/all_cuda-129_arch-x86_64.yaml +++ b/conda/environments/all_cuda-129_arch-x86_64.yaml @@ -61,6 +61,7 @@ dependencies: - requests - rmm==26.2.*,>=0.0.0a0 - scikit-build-core>=0.11.0 +- scipy>=1.13.0 - sphinx - sphinx-copybutton - sphinx-design diff --git a/conda/environments/all_cuda-131_arch-aarch64.yaml b/conda/environments/all_cuda-131_arch-aarch64.yaml index fb23f887a..017c779ee 100644 --- a/conda/environments/all_cuda-131_arch-aarch64.yaml +++ b/conda/environments/all_cuda-131_arch-aarch64.yaml @@ -61,6 +61,7 @@ dependencies: - requests - rmm==26.2.*,>=0.0.0a0 - scikit-build-core>=0.11.0 +- scipy>=1.13.0 - sphinx - sphinx-copybutton - sphinx-design diff --git a/conda/environments/all_cuda-131_arch-x86_64.yaml b/conda/environments/all_cuda-131_arch-x86_64.yaml index 501729acd..71f480f36 100644 --- a/conda/environments/all_cuda-131_arch-x86_64.yaml +++ b/conda/environments/all_cuda-131_arch-x86_64.yaml @@ -61,6 +61,7 @@ dependencies: - requests - rmm==26.2.*,>=0.0.0a0 - scikit-build-core>=0.11.0 +- scipy>=1.13.0 - sphinx - sphinx-copybutton - sphinx-design diff --git a/conda/recipes/cuopt/recipe.yaml b/conda/recipes/cuopt/recipe.yaml index 3408c6715..7c9ecf4f7 100644 --- a/conda/recipes/cuopt/recipe.yaml +++ b/conda/recipes/cuopt/recipe.yaml @@ -82,6 +82,7 @@ requirements: - pylibraft =${{ minor_version }} - python - rmm =${{ minor_version }} + - scipy # Needed by Numba for CUDA support - cuda-nvcc-impl # TODO: Add nvjitlink here diff --git a/dependencies.yaml b/dependencies.yaml index 7dc6b9490..b3da7c656 100644 --- a/dependencies.yaml +++ b/dependencies.yaml @@ -350,6 +350,7 @@ dependencies: - numba>=0.60.0 - &pandas pandas>=2.0 - &pyyaml pyyaml>=6.0.0 + - scipy>=1.13.0 - output_types: requirements packages: # pip recognizes the index as a global option for the requirements.txt file diff --git a/python/cuopt/cuopt/linear_programming/problem.py b/python/cuopt/cuopt/linear_programming/problem.py index 4259753dc..f5b9ba2b2 100644 --- a/python/cuopt/cuopt/linear_programming/problem.py +++ b/python/cuopt/cuopt/linear_programming/problem.py @@ -1,4 +1,4 @@ -# SPDX-FileCopyrightText: Copyright (c) 2025 NVIDIA CORPORATION & AFFILIATES. All rights reserved. +# SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. # SPDX-License-Identifier: Apache-2.0 @@ -7,6 +7,7 @@ import cuopt_mps_parser import numpy as np +from scipy.sparse import coo_matrix import cuopt.linear_programming.data_model as data_model import cuopt.linear_programming.solver as solver @@ -326,7 +327,13 @@ class QuadraticExpression: """ def __init__( - self, qvars1, qvars2, qcoefficients, vars, coefficients, constant + self, + qvars1, + qvars2, + qcoefficients, + vars=[], + coefficients=[], + constant=0.0, ): self.qvars1 = qvars1 self.qvars2 = qvars2 @@ -681,6 +688,204 @@ def __eq__(self, other): raise Exception("Quadratic constraints not supported") +class MQuadraticExpression: + """ + MQuadraticExpressions contain quadratic terms, linear terms, and a constant. + MQuadraticExpressions can be used to create quadratic objectives in + the Problem. + + Parameters + ---------- + qmatrix : List[List[float]] + 2D List or np.array containing quadratic coefficient matrix terms. + qvars : List[Variable] + List of variables for quadratic matrix. + Should be in the same order of variables added. + vars : List[Variable] + List of Variables for linear terms. + coefficients : List[float] + List of coefficients for linear terms. + constant : float + Constant of the quadratic expression. + + Examples + -------- + >>> x = problem.addVariable() + >>> y = problem.addVariable() + >>> # Create x^2 + 2*x*y + 3*x + 4 + >>> quad_expr = MQuadraticExpression( + ... [[1.0, 2.0], [0.0, 0.0]], [x, y], + ... [x], [3.0], 4.0 + ... ) + """ + + def __init__( + self, qmatrix, qvars=[], vars=[], coefficients=[], constant=0.0 + ): + self.qmatrix = qmatrix + self.qvars = qvars + self.vars = vars + self.coefficients = coefficients + self.constant = constant + + def shape(self): + return np.shape(self.qmatrix) + + def getValue(self): + """ + Returns the value of the expression computed with the + current solution. + """ + # TODO: improve efficiency + value = 0.0 + for i, var in enumerate(self.vars): + value += var.Value * self.coefficients[i] + for i, var1 in enumerate(self.qvars): + for j, var2 in enumerate(self.qvars): + value += var1.Value * var2.Value * float(self.qmatrix[i, j]) + return value + self.constant + + def __size__(self): + return np.size(self.qmatrix) + + def __iadd__(self, other): + # Compute quad_expr += lin_expr + match other: + case int() | float(): + # Update just the constant value + self.constant += float(other) + return self + case Variable(): + # Append just a variable with coefficient 1.0 + self.vars.append(other) + self.coefficients.append(1.0) + return self + case LinearExpression(): + # Append all linear variables, coefficients and constants + self.vars.extend(other.vars) + self.coefficients.extend(other.coefficients) + self.constant += other.constant + return self + + def __add__(self, other): + # Compute quad_expr1 = quar_expr2 + lin_expr + match other: + case int() | float(): + # Update just the constant value + return MQuadraticExpression( + self.qmatrix, + self.qvars, + self.vars, + self.coefficients, + self.constant + float(other), + ) + case Variable(): + # Append just a variable with coefficient 1.0 + vars = self.vars + [other] + coeffs = self.coefficients + [1.0] + return MQuadraticExpression( + self.qmatrix, + self.qvars, + vars, + coeffs, + self.constant, + ) + case LinearExpression(): + # Append all linear variables, coefficients and constants + vars = self.vars + other.vars + coeffs = self.coefficients + other.coefficients + constant = self.constant + other.constant + return MQuadraticExpression( + self.qmatrix, + self.qvars, + vars, + coeffs, + constant, + ) + + def __isub__(self, other): + # Compute quad_expr -= lin_expr + match other: + case int() | float(): + # Update just the constant value + self.constant -= float(other) + return self + case Variable(): + # Append just a variable with coefficient -1.0 + self.vars.append(other) + self.coefficients.append(-1.0) + return self + case LinearExpression(): + # Append all linear variables, coefficients and constants + self.vars.extend(other.vars) + for coeff in other.coefficients: + self.coefficients.append(-coeff) + self.constant -= other.constant + return self + + def __sub__(self, other): + # Compute quad_expr2 = quad_expr1 - lin_expr + match other: + case int() | float(): + # Update just the constant value + return MQuadraticExpression( + self.qmatrix, + self.qvars, + self.vars, + self.coefficients, + self.constant - float(other), + ) + case Variable(): + # Append just a variable with coefficient -1.0 + vars = self.vars + [other] + coeffs = self.coefficients + [-1.0] + return MQuadraticExpression( + self.qmatrix, + self.qvars, + vars, + coeffs, + self.constant, + ) + case LinearExpression(): + # Append all linear variables, coefficients and constants + vars = self.vars + other.vars + coeffs = [] + for i in self.coefficients: + coeffs.append(i) + for i in other.coefficients: + coeffs.append(-1.0 * i) + constant = self.constant - other.constant + return MQuadraticExpression( + self.qmatrix, + self.qvars, + vars, + coeffs, + constant, + ) + + ## TODO: Add matrix multiplication + # def __matmul__(self, qcols): + # if not self.qcols: + # self.qcols = qcols + # else: + # raise Exception("") + + # def __rmatmul__(self, qcols): + # if not self.qrows: + # self.qrows = qrows + # else: + # raise Exception("") + + def __le__(self, other): + raise Exception("Quadratic constraints not supported") + + def __ge__(self, other): + raise Exception("Quadratic constraints not supported") + + def __eq__(self, other): + raise Exception("Quadratic constraints not supported") + + class LinearExpression: """ LinearExpressions contain a set of variables, the coefficients @@ -1140,7 +1345,6 @@ def __init__(self, model_name=""): self.ObjSense = MINIMIZE self.ObjConstant = 0.0 self.Status = -1 - self.ObjValue = float("nan") self.warmstart_data = None self.model = None @@ -1149,7 +1353,7 @@ def __init__(self, model_name=""): self.row_sense = None self.constraint_csr_matrix = None self.objective_qcsr_matrix = None - self.objective_qcoo_matrix = [], [], [] + self.objective_qcoo_matrix = None self.lower_bound = None self.upper_bound = None self.var_type = None @@ -1277,9 +1481,9 @@ def _to_data_model(self): dm.set_objective_offset(self.ObjConstant) if self.getQcsr(): dm.set_quadratic_objective_matrix( - np.array(self.objective_qcsr_matrix["values"]), - np.array(self.objective_qcsr_matrix["column_indices"]), - np.array(self.objective_qcsr_matrix["row_pointers"]), + self.objective_qcsr_matrix["values"], + self.objective_qcsr_matrix["column_indices"], + self.objective_qcsr_matrix["row_pointers"], ) dm.set_variable_lower_bounds(self.lower_bound) dm.set_variable_upper_bounds(self.upper_bound) @@ -1310,9 +1514,8 @@ def reset_solved_values(self): self.model = None self.constraint_csr_matrix = None - self.objective_qcoo_matrix = [], [], [] + self.objective_qcoo_matrix = None self.objective_qcsr_matrix = None - self.ObjValue = float("nan") self.warmstart_data = None self.solved = False @@ -1480,11 +1683,26 @@ def setObjective(self, expr, sense=MINIMIZE): sum_coeff ) self.ObjConstant = expr.constant - self.objective_qcoo_matrix = ( - expr.qvars1, - expr.qvars2, - expr.qcoefficients, + qrows = [var.getIndex() for var in expr.qvars1] + qcols = [var.getIndex() for var in expr.qvars2] + self.objective_qcoo_matrix = coo_matrix( + ( + np.array(expr.qcoefficients), + (np.array(qrows), np.array(qcols)), + ), + shape=(self.NumVariables, self.NumVariables), ) + case MQuadraticExpression(): + for var in self.vars: + var.setObjectiveCoefficient(0.0) + for var, coeff in zip(expr.vars, expr.coefficients): + c_val = self.vars[var.getIndex()].getObjectiveCoefficient() + sum_coeff = coeff + c_val + self.vars[var.getIndex()].setObjectiveCoefficient( + sum_coeff + ) + self.ObjConstant = expr.constant + self.objective_qcoo_matrix = coo_matrix(expr.qmatrix) case _: raise ValueError( "Objective must be a LinearExpression or a constant" @@ -1643,16 +1861,21 @@ def Obj(self): coeffs = [] for var in self.vars: coeffs.append(var.getObjectiveCoefficient()) - if not self.objective_qcoo_matrix: + if self.objective_qcoo_matrix is None: return LinearExpression(self.vars, coeffs, self.ObjConstant) else: - return QuadraticExpression( - *self.objective_qcoo_matrix, + return MQuadraticExpression( + self.objective_qcoo_matrix.todense(), + self.vars, self.vars, coeffs, self.ObjConstant, ) + @property + def ObjValue(self): + return self.Obj.getValue() + def getCSR(self): """ Computes and returns the CSR representation of the @@ -1673,29 +1896,15 @@ def getCSR(self): def getQcsr(self): if self.objective_qcsr_matrix is not None: return self.dict_to_object(self.objective_qcsr_matrix) - vars1, vars2, coeffs = self.objective_qcoo_matrix - if not vars1: + if self.objective_qcoo_matrix is None: return None - Qdict = {} - Qcsr_dict = {"row_pointers": [0], "column_indices": [], "values": []} - for i, var1 in enumerate(vars1): - if var1.index not in Qdict: - Qdict[var1.index] = {} - row_dict = Qdict[var1.index] - var2 = vars2[i] - coeff = coeffs[i] - row_dict[var2.index] = ( - row_dict[var2.index] + coeff - if var2.index in row_dict - else coeff - ) - for i in range(0, self.NumVariables): - if i in Qdict: - Qcsr_dict["column_indices"].extend(list(Qdict[i].keys())) - Qcsr_dict["values"].extend(list(Qdict[i].values())) - Qcsr_dict["row_pointers"].append(len(Qcsr_dict["column_indices"])) - self.objective_qcsr_matrix = Qcsr_dict - return self.dict_to_object(Qcsr_dict) + qcsr_matrix = self.objective_qcoo_matrix.tocsr() + self.objective_qcsr_matrix = { + "row_pointers": qcsr_matrix.indptr, + "column_indices": qcsr_matrix.indices, + "values": qcsr_matrix.data, + } + return self.dict_to_object(self.objective_qcsr_matrix) def relax(self): """ @@ -1725,7 +1934,6 @@ def populate_solution(self, solution): else: IsMIP = True self.SolutionStats = self.dict_to_object(solution.get_milp_stats()) - primal_sol = solution.get_primal_solution() reduced_cost = solution.get_reduced_cost() if len(primal_sol) > 0: @@ -1744,7 +1952,6 @@ def populate_solution(self, solution): if dual_sol is not None and len(dual_sol) > 0: constr.DualValue = dual_sol[i] constr.Slack = constr.compute_slack() - self.ObjValue = self.Obj.getValue() self.solved = True def solve(self, settings=solver_settings.SolverSettings()): @@ -1765,9 +1972,7 @@ def solve(self, settings=solver_settings.SolverSettings()): """ if self.model is None: self._to_data_model() - # Call Solver solution = solver.Solve(self.model, settings) - # Post Solve self.populate_solution(solution) diff --git a/python/cuopt/cuopt/tests/linear_programming/test_python_API.py b/python/cuopt/cuopt/tests/linear_programming/test_python_API.py index 06ed93703..b01d43611 100644 --- a/python/cuopt/cuopt/tests/linear_programming/test_python_API.py +++ b/python/cuopt/cuopt/tests/linear_programming/test_python_API.py @@ -1,4 +1,4 @@ -# SPDX-FileCopyrightText: Copyright (c) 2025 NVIDIA CORPORATION & AFFILIATES. All rights reserved. +# SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. # SPDX-License-Identifier: Apache-2.0 import math @@ -20,6 +20,7 @@ Problem, VType, sense, + MQuadraticExpression, ) from cuopt.linear_programming.solver.solver_parameters import ( CUOPT_AUGMENTED, @@ -754,9 +755,9 @@ def test_quadratic_expression_and_matrix(): exp_col_inds = [0, 1, 2, 0, 1, 2] exp_vals = [21, -7, 7, 3, -1, 1] - assert Qcsr.row_pointers == exp_row_ptrs - assert Qcsr.column_indices == exp_col_inds - assert Qcsr.values == exp_vals + assert list(Qcsr.row_pointers) == exp_row_ptrs + assert list(Qcsr.column_indices) == exp_col_inds + assert list(Qcsr.values) == exp_vals def test_quadratic_objective_1(): @@ -813,3 +814,35 @@ def test_quadratic_objective_2(): assert x2.getValue() == pytest.approx(0.0000000, abs=0.000001) assert x3.getValue() == pytest.approx(0.1092896, abs=1e-3) assert problem.ObjValue == pytest.approx(-0.284153, abs=1e-3) + + +def test_quadratic_matrix(): + # Minimize 4 x1^2 + 2 x2^2 + 3 x3^2 + 1.5 x1 x3 - 2 x1 + 0.5 x2 - x3 + 4 + # subject to x1 + 2*x2 + x3 <= 3 + # x1 >= 0 + # x2 >= 0 + # x3 >= 0 + + problem = Problem() + x1 = problem.addVariable(lb=0, name="x") + x2 = problem.addVariable(lb=0, name="y") + x3 = problem.addVariable(lb=0, name="z") + + problem.addConstraint(x1 + 2 * x2 + x3 <= 3) + + Q = [[4, 0, 1.5], [0, 2, 0], [0, 0, 3]] + quad_expr = MQuadraticExpression(Q) + quad_expr1 = quad_expr + 4 # Quad_matrix add constant + quad_expr2 = quad_expr1 - x3 # Quad_matrix sub variable + quad_expr2 -= 2 * x1 # Quad_matrix isub lin_expr + quad_expr2 += 0.5 * x2 # Quad_matrix iadd lin_expr + + problem.setObjective(quad_expr2) + + problem.solve() + + assert problem.Status.name == "Optimal" + assert x1.getValue() == pytest.approx(0.2295081, abs=1e-3) + assert x2.getValue() == pytest.approx(0.0000000, abs=0.000001) + assert x3.getValue() == pytest.approx(0.1092896, abs=1e-3) + assert problem.ObjValue == pytest.approx(3.715847, abs=1e-3) diff --git a/python/cuopt/pyproject.toml b/python/cuopt/pyproject.toml index 05ca7ffa8..0e8eedc64 100644 --- a/python/cuopt/pyproject.toml +++ b/python/cuopt/pyproject.toml @@ -32,6 +32,7 @@ dependencies = [ "pyyaml>=6.0.0", "rapids-logger==0.2.*,>=0.0.0a0", "rmm==26.2.*,>=0.0.0a0", + "scipy>=1.13.0", ] # This list was generated by `rapids-dependency-file-generator`. To make changes, edit ../../dependencies.yaml and run `rapids-dependency-file-generator`. classifiers = [ "Intended Audience :: Developers",