diff --git a/python/cuopt/cuopt/linear_programming/data_model/data_model_wrapper.pyx b/python/cuopt/cuopt/linear_programming/data_model/data_model_wrapper.pyx index ec8bdf3730..b298cf57a5 100644 --- a/python/cuopt/cuopt/linear_programming/data_model/data_model_wrapper.pyx +++ b/python/cuopt/cuopt/linear_programming/data_model/data_model_wrapper.pyx @@ -1,4 +1,4 @@ -# SPDX-FileCopyrightText: Copyright (c) 2023-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. # noqa +# SPDX-FileCopyrightText: Copyright (c) 2023-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. # SPDX-License-Identifier: Apache-2.0 @@ -292,6 +292,9 @@ cdef class DataModel: return self.problem_name def set_data_model_view(self): + # Rebind from scratch so optional fields that were previously set but + # are now empty (e.g. Q after QP -> LP) do not stick in the C++ view. + self.c_data_model_view.reset(new data_model_view_t[int, double]()) cdef data_model_view_t[int, double]* c_data_model_view = ( self.c_data_model_view.get() ) diff --git a/python/cuopt/cuopt/linear_programming/problem.py b/python/cuopt/cuopt/linear_programming/problem.py index 49885acf52..2c495af53c 100644 --- a/python/cuopt/cuopt/linear_programming/problem.py +++ b/python/cuopt/cuopt/linear_programming/problem.py @@ -5,9 +5,10 @@ import copy import os from enum import Enum +from itertools import chain import numpy as np -from scipy.sparse import coo_matrix +from scipy.sparse import coo_matrix, csr_matrix import cuopt.linear_programming.data_model as data_model from cuopt.linear_programming import ParseMps, Read @@ -111,6 +112,17 @@ class Variable: (unset). Only used for a MIP problem. """ + # Input attrs: direct writes mark the matching Problem._stale key. + _INPUT_ATTRIBUTE = { + "LB": "variable", + "UB": "variable", + "Obj": "objective", + "VariableType": "variable", + "VariableName": "variable", + "MIPStart": "variable", + } + _OUTPUT_ATTRIBUTES = frozenset({"Value", "ReducedCost"}) + def __init__( self, lb=0.0, @@ -129,6 +141,18 @@ def __init__( self.VariableName = vname self.MIPStart = float("nan") + def __setattr__(self, name, value): + object.__setattr__(self, name, value) + # Solution fields are written in hot loops; skip tracking lookups. + if name in self._OUTPUT_ATTRIBUTES: + return + problem = self.__dict__.get("_problem") + stale_key = self._INPUT_ATTRIBUTE.get(name) + if problem is not None and stale_key is not None: + if problem.solved: + problem._reset_solved_values() + problem._mark_stale(stale_key) + def getIndex(self): """ Get the index position of the variable in the problem. @@ -1311,11 +1335,20 @@ class Constraint: is_quadratic : bool True when the row is exported as a QCMATRIX quadratic constraint. Slack : float - Computed LHS - RHS with current solution. + Classical LP slack/surplus for the current solution: ``rhs - lhs`` + for ``<=``, ``lhs - rhs`` for ``>=`` (non-negative if feasible). For + ``==``, the residual ``rhs - lhs`` (near zero if feasible). DualValue : float Constraint dual value in the current solution. """ + _INPUT_ATTRIBUTE = { + "RHS": "rhs", + "Sense": "structure", + "ConstraintName": "structure", + } + _OUTPUT_ATTRIBUTES = frozenset({"Slack", "DualValue"}) + def __init__(self, expr, sense, rhs, name=""): self.index = -1 self.Sense = sense @@ -1358,6 +1391,21 @@ def __init__(self, expr, sense, rhs, name=""): ) self.RHS = rhs - expr.getConstant() + def __setattr__(self, name, value): + object.__setattr__(self, name, value) + # Solution fields are written in hot loops; skip tracking lookups. + if name in self._OUTPUT_ATTRIBUTES: + return + problem = self.__dict__.get("_problem") + stale_key = self._INPUT_ATTRIBUTE.get(name) + if problem is not None and stale_key is not None: + if name == "RHS" and self.is_quadratic: + object.__setattr__(self, "rhs_value", value) + stale_key = "structure" + if problem.solved: + problem._reset_solved_values() + problem._mark_stale(stale_key) + def __len__(self): return len(self.vindex_coeff_dict) @@ -1387,16 +1435,6 @@ def getCoefficient(self, var): v_idx = var.index return self.vindex_coeff_dict[v_idx] - def compute_slack(self): - # Computes the constraint Slack in the current solution. - index_to_var = {var.index: var for var in self.vars} - lhs = sum( - index_to_var[v_idx].Value * coeff - for v_idx, coeff in self.vindex_coeff_dict.items() - ) - - return self.RHS - lhs - class Problem: """ @@ -1466,6 +1504,32 @@ def __init__(self, model_name=""): self.lower_bound = None self.upper_bound = None self.var_type = None + self._constraint_index_to_csr_row = None + # Value/structure dirtiness for warm DataModel sync. + # True = must refresh that category before the next solve. + self._stale = { + "structure": True, + "variable": True, + "objective": True, + "rhs": True, + "A_values": True, + } + + def _mark_stale(self, *keys): + for key in keys: + if key not in self._stale: + raise KeyError(f"Unknown stale key: {key}") + self._stale[key] = True + + def _clear_stale(self, *keys): + if not keys: + for key in self._stale: + self._stale[key] = False + return + for key in keys: + if key not in self._stale: + raise KeyError(f"Unknown stale key: {key}") + self._stale[key] = False class dict_to_object: def __init__(self, mdict): @@ -1523,48 +1587,79 @@ def _from_data_model(self, dm): else: raise Exception("Couldn't initialize constraints") - def _to_data_model(self): - dm = data_model.DataModel() - - # iterate through the constraints and construct the constraint matrix - n = len(self.vars) - self.rhs = [] - self.row_sense = [] - self.row_names = [] - - if self.constraint_csr_matrix is None: - csr_dict = { - "row_pointers": [0], - "column_indices": [], - "values": [], - } - for constr in self.constrs: - if constr.is_quadratic: - continue - csr_dict["column_indices"].extend( - list(constr.vindex_coeff_dict.keys()) - ) - csr_dict["values"].extend( - list(constr.vindex_coeff_dict.values()) - ) - csr_dict["row_pointers"].append( - len(csr_dict["column_indices"]) - ) - self.rhs.append(constr.RHS) - self.row_sense.append(constr.Sense) - constr_name = constr.ConstraintName - if constr_name == "": - constr_name = "R" + str(constr.index) - self.row_names.append(constr_name) - self.constraint_csr_matrix = csr_dict - + # Quadratic (QCMATRIX) rows are kept out of the linear CSR by the + # reader, each in its own bundle. Rebuild them as quadratic + # Constraints so the Python problem is the whole model that was read. + for qc in dm.get_quadratic_constraints(): + expr = QuadraticExpression( + qvars1=[vars[i] for i in qc["rows"]], + qvars2=[vars[j] for j in qc["cols"]], + qcoefficients=qc["vals"], + vars=[vars[j] for j in qc["linear_indices"]], + coefficients=qc["linear_values"], + ) + rhs = qc["rhs_value"] + row = ( + expr <= rhs if qc["constraint_row_type"] == LE else expr >= rhs + ) + self.addConstraint(row, name=qc["constraint_row_name"]) + + # setObjective(linear) leaves objective_qmatrix as None. Copy Q from + # the source DataModel so later _to_data_model / writeMPS / solve keep + # the quadratic objective. + Q_values = dm.get_quadratic_objective_values() + Q_indices = dm.get_quadratic_objective_indices() + Q_offsets = dm.get_quadratic_objective_offsets() + if len(Q_values) > 0: + self.objective_qmatrix = csr_matrix( + (Q_values, Q_indices, Q_offsets), + shape=(num_vars, num_vars), + ) else: - for constr in self.constrs: - if constr.is_quadratic: - continue - self.rhs.append(constr.RHS) - self.row_sense.append(constr.Sense) + self.objective_qmatrix = None + + # Adopt the source CSR and value arrays instead of leaving the caches + # empty. Without them the first solve/writeMPS falls back to a full + # rebuild from Python, which costs an O(nnz) pass and drops whatever + # the reader set but Python does not model. + self._rebuild_row_caches() + self._rebuild_variable_caches() + self.constraint_csr_matrix = { + "row_pointers": np.array(offsets, dtype=np.int32), + "column_indices": np.array(indices, dtype=np.int32), + "values": np.array(values, dtype=np.float64), + } + def _rebuild_row_caches(self): + """Refresh the per-row caches and return the linear constraints. + + Quadratic rows are excluded: they are not part of the linear CSR and + carry their own RHS, so every cache here is indexed by CSR row. + """ + linear_constrs = [ + constr for constr in self.constrs if not constr.is_quadratic + ] + m = len(linear_constrs) + self._constraint_index_to_csr_row = { + constr.index: row for row, constr in enumerate(linear_constrs) + } + self.rhs = np.fromiter( + (constr.RHS for constr in linear_constrs), + dtype=np.float64, + count=m, + ) + self.row_sense = np.asarray( + [constr.Sense for constr in linear_constrs], dtype="S1" + ) + self.row_names = [ + constr.ConstraintName or "R" + str(constr.index) + for constr in linear_constrs + ] + return linear_constrs + + def _rebuild_variable_caches(self): + """Refresh the per-variable caches from the Variable objects.""" + n = len(self.vars) self.objective = np.zeros(n) self.lower_bound, self.upper_bound = np.zeros(n), np.zeros(n) self.var_type = np.empty(n, dtype="S1") @@ -1582,16 +1677,56 @@ def _to_data_model(self): self.var_names.append(var_name) self.mip_start[j] = self.vars[j].MIPStart + def _to_data_model(self): + dm = data_model.DataModel() + + # A full DataModel rebuild always regenerates CSR from the constraint + # dictionaries. Warm value-only paths use _refresh_data_model_values(). + linear_constrs = self._rebuild_row_caches() + m = len(linear_constrs) + + row_sizes = np.fromiter( + (len(constr.vindex_coeff_dict) for constr in linear_constrs), + dtype=np.int32, + count=m, + ) + row_pointers = np.empty(m + 1, dtype=np.int32) + row_pointers[0] = 0 + np.cumsum(row_sizes, out=row_pointers[1:]) + nnz = int(row_pointers[-1]) + + column_indices = np.fromiter( + chain.from_iterable( + constr.vindex_coeff_dict.keys() for constr in linear_constrs + ), + dtype=np.int32, + count=nnz, + ) + values = np.fromiter( + chain.from_iterable( + constr.vindex_coeff_dict.values() for constr in linear_constrs + ), + dtype=np.float64, + count=nnz, + ) + self.constraint_csr_matrix = { + "row_pointers": row_pointers, + "column_indices": column_indices, + "values": values, + } + + self._rebuild_variable_caches() + # Initialize datamodel dm.set_csr_constraint_matrix( - np.array(self.constraint_csr_matrix["values"]), - np.array(self.constraint_csr_matrix["column_indices"]), - np.array(self.constraint_csr_matrix["row_pointers"]), + self.constraint_csr_matrix["values"], + self.constraint_csr_matrix["column_indices"], + self.constraint_csr_matrix["row_pointers"], ) if self.ObjSense == -1: dm.set_maximize(True) - dm.set_constraint_bounds(np.array(self.rhs)) - dm.set_row_types(np.array(self.row_sense, dtype="S1")) + dm.set_constraint_bounds(self.rhs) + dm.set_row_types(self.row_sense) dm.set_objective_coefficients(self.objective) dm.set_objective_offset(self.ObjConstant) if self.objective_qmatrix is not None: @@ -1628,6 +1763,106 @@ def _to_data_model(self): dm.set_initial_primal_solution(self.mip_start) self.model = dm + self._clear_stale() + + def _invalidate_problem_cache(self): + """Drop DataModel/CSR caches and mark all categories stale. + + Used when variable/constraint topology changes. Does not clear + objective_qmatrix (Q lifetime is independent of A structure). + """ + self.model = None + self.constraint_csr_matrix = None + self._constraint_index_to_csr_row = None + self._mark_stale( + "structure", "variable", "objective", "rhs", "A_values" + ) + + def _refresh_data_model_values(self): + """Patch existing DataModel when sparsity structure is unchanged.""" + n = len(self.vars) + if ( + self.model is None + or self.constraint_csr_matrix is None + or self._stale["structure"] + ): + self._to_data_model() + return + + # Nothing changed since last sync. + if not any(self._stale.values()): + return + + if self._stale["objective"]: + for j in range(n): + self.objective[j] = self.vars[j].getObjectiveCoefficient() + self.model.set_objective_coefficients(self.objective) + self.model.set_objective_offset(self.ObjConstant) + self.model.set_maximize(self.ObjSense == -1) + # Q is independent of A topology; replace or clear in place. + # Empty arrays clear Python Q_*; set_data_model_view rebuilds the + # C++ view each Solve, so a prior Q does not stick (QP -> LP). + if self.objective_qmatrix is not None: + self.model.set_quadratic_objective_matrix( + self.objective_qmatrix.data, + self.objective_qmatrix.indices, + self.objective_qmatrix.indptr, + ) + else: + self.model.set_quadratic_objective_matrix( + np.array([], dtype=np.float64), + np.array([], dtype=np.int32), + np.array([], dtype=np.int32), + ) + self._stale["objective"] = False + + if self._stale["variable"]: + # Rebuild non-objective variable arrays from Variable objects. + self.var_names = [] + for j in range(n): + var = self.vars[j] + self.var_type[j] = var.getVariableType() + self.lower_bound[j] = var.getLowerBound() + self.upper_bound[j] = var.getUpperBound() + var_name = var.VariableName + if var_name == "": + var_name = "C" + str(var.index) + self.var_names.append(var_name) + self.mip_start[j] = var.MIPStart + self.model.set_variable_lower_bounds(self.lower_bound) + self.model.set_variable_upper_bounds(self.upper_bound) + self.model.set_variable_types(self.var_type) + self.model.set_variable_names(self.var_names) + if self.mip_start.size > 0 and not np.all( + np.isnan(self.mip_start) + ): + self.model.set_initial_primal_solution(self.mip_start) + else: + # Empty so the next Solve skips add_initial_mip_solution / + # set_initial_pdlp_primal_solution (clears a prior warm start). + self.model.set_initial_primal_solution(np.array([])) + self._stale["variable"] = False + + if self._stale["rhs"]: + # self.rhs is a float64 ndarray after _to_data_model(). + linear_row = 0 + for constr in self.constrs: + if constr.is_quadratic: + continue + self.rhs[linear_row] = constr.RHS + linear_row += 1 + self.model.set_constraint_bounds(self.rhs) + self._stale["rhs"] = False + + if self._stale["A_values"]: + # Existing nonzero values changed but CSR topology is stable. + # Reuse indices/offsets and synchronize the updated values once. + self.model.set_csr_constraint_matrix( + self.constraint_csr_matrix["values"], + self.constraint_csr_matrix["column_indices"], + self.constraint_csr_matrix["row_pointers"], + ) + self._stale["A_values"] = False def update(self): """ @@ -1635,24 +1870,34 @@ def update(self): existing Variables, Constraints or Objective has been modified. """ - self.reset_solved_values() + self._reset_solved_values() + # Direct attribute edits (e.g. x.Obj / c.RHS) are opaque; mark all + # value categories stale. Drop CSR because direct coefficient edits + # cannot be mapped back to cached A_values reliably. + self.constraint_csr_matrix = None + self._mark_stale("variable", "objective", "rhs", "A_values") - def reset_solved_values(self): - # Resets all post solve values - for var in self.vars: - var.Value = float("nan") - var.ReducedCost = float("nan") + def _reset_solved_values( + self, *, invalidate_structure=False, invalidate_solution=True + ): + # Resets post-solve values and/or drops structure caches. + if invalidate_solution: + for var in self.vars: + var.Value = float("nan") + var.ReducedCost = float("nan") - for constr in self.constrs: - constr.Slack = float("nan") - constr.DualValue = float("nan") + for constr in self.constrs: + constr.Slack = float("nan") + constr.DualValue = float("nan") - self.model = None - self.constraint_csr_matrix = None - self.objective_qmatrix = None self.warmstart_data = None self.solved = False + if invalidate_structure: + # Constraint topology only; the quadratic objective is unrelated + # and must survive. + self._invalidate_problem_cache() + def addVariable( self, lb=0.0, ub=float("inf"), obj=0.0, vtype=CONTINUOUS, name="" ): @@ -1684,11 +1929,16 @@ def addVariable( name="Var1") """ if self.solved: - self.reset_solved_values() # Reset all solved values + self._reset_solved_values() + self._invalidate_problem_cache() n = len(self.vars) var = Variable(lb, ub, obj, vtype, name) var.index = n + var._problem = self self.vars.append(var) + if self.objective_qmatrix is not None: + # Q is sized by variable count; the new column/row is all zeros. + self.objective_qmatrix.resize((n + 1, n + 1)) return var def addConstraint(self, constr, name=""): @@ -1715,12 +1965,14 @@ def addConstraint(self, constr, name=""): >>> problem.addConstraint(-x*x + y*y <= 0, name="soc") """ if self.solved: - self.reset_solved_values() # Reset all solved values + self._reset_solved_values() + self._invalidate_problem_cache() n = len(self.constrs) match constr: case Constraint(): constr.index = n constr.ConstraintName = name + constr._problem = self self.constrs.append(constr) case _: raise ValueError("addConstraint requires a Constraint object") @@ -1750,23 +2002,70 @@ def updateConstraint(self, constr, coeffs=None, rhs=None): >>> c2 = problem.addConstraint(x + y <= 5, name="c2") >>> problem.updateConstraint(c1, coeffs=[(x, 1)], rhs=10) """ - self.reset_solved_values() if coeffs is None: coeffs = [] - if isinstance(constr, Constraint): - if constr.is_quadratic: - raise ValueError( - "updateConstraint applies to linear constraints only" - ) - if isinstance(coeffs, dict): - coeffs = coeffs.items() - for var, coeff in coeffs: - idx = var.index - constr.vindex_coeff_dict[idx] = coeff - if rhs is not None: - constr.RHS = rhs - else: + if not isinstance(constr, Constraint): raise ValueError("Object to update must be a Constraint") + if ( + constr.index < 0 + or constr.index >= len(self.constrs) + or self.constrs[constr.index] is not constr + ): + raise ValueError("Constraint does not belong to this Problem") + if constr.is_quadratic: + raise ValueError( + "updateConstraint applies to linear constraints only" + ) + if isinstance(coeffs, dict): + coeffs = list(coeffs.items()) + else: + coeffs = list(coeffs) + + has_new_nonzero = False + if coeffs: + for var, coeff in coeffs: + if var.index not in constr.vindex_coeff_dict: + has_new_nonzero = True + constr.vindex_coeff_dict[var.index] = coeff + + if not has_new_nonzero: + if ( + self.constraint_csr_matrix is not None + and not self._stale["structure"] + ): + row = self._constraint_index_to_csr_row[constr.index] + row_pointers = self.constraint_csr_matrix["row_pointers"] + column_indices = self.constraint_csr_matrix[ + "column_indices" + ] + values = self.constraint_csr_matrix["values"] + start = int(row_pointers[row]) + end = int(row_pointers[row + 1]) + + for var, coeff in coeffs: + position = start + int( + np.flatnonzero( + column_indices[start:end] == var.index + )[0] + ) + values[position] = coeff + + # Defer DataModel sync until solve() so multiple edits batch. + self._mark_stale("A_values") + + if rhs is not None: + constr.RHS = rhs + self._mark_stale("rhs") + + if has_new_nonzero: + # Topology changed: always drop structure caches. Skip the O(n) + # solution wipe when already dirty from a prior mutation. + self._reset_solved_values( + invalidate_structure=True, + invalidate_solution=self.solved, + ) + elif self.solved: + self._reset_solved_values() def setObjective(self, expr, sense=MINIMIZE): """ @@ -1793,8 +2092,10 @@ def setObjective(self, expr, sense=MINIMIZE): >>> problem.setObjective(x + y, sense=MAXIMIZE) """ if self.solved: - self.reset_solved_values() # Reset all solved values + self._reset_solved_values() # Reset all solved values self.ObjSense = sense + self._mark_stale("objective") + is_quadratic = False match expr: case int() | float(): for var in self.vars: @@ -1816,6 +2117,7 @@ def setObjective(self, expr, sense=MINIMIZE): ) self.ObjConstant = expr.getConstant() case QuadraticExpression(): + is_quadratic = True for var in self.vars: var.setObjectiveCoefficient(0.0) for var, coeff in zip(expr.vars, expr.coefficients): @@ -1841,6 +2143,10 @@ def setObjective(self, expr, sense=MINIMIZE): raise ValueError( "Objective must be a Variable, Expression or a constant" ) + if not is_quadratic: + # QP -> LP: drop Python Q; warm refresh pushes empty Q and + # set_data_model_view rebuilds the C++ view so old Q is gone. + self.objective_qmatrix = None def updateObjective(self, coeffs=None, constant=None, sense=None): """ @@ -1866,7 +2172,8 @@ def updateObjective(self, coeffs=None, constant=None, sense=None): >>> problem.updateObjective(coeffs=[(x1, 1.0), (x2, 3.0)], constant=5, sense=MINIMIZE) """ - self.reset_solved_values() + if self.solved: + self._reset_solved_values() if coeffs is None: coeffs = [] if isinstance(coeffs, dict): @@ -1877,6 +2184,7 @@ def updateObjective(self, coeffs=None, constant=None, sense=None): self.ObjConstant = constant if sense is not None: self.ObjSense = sense + self._mark_stale("objective") def getIncumbentValues(self, solution, vars): """ @@ -2002,6 +2310,9 @@ def read(cls, file_path, fixed_mps_format=False): data_model = Read(file_path, fixed_mps_format) problem._from_data_model(data_model) problem.model = data_model + # Python objects match the attached file DataModel; avoid a no-op + # structure rebuild that would rewrite the model from Python. + problem._clear_stale() return problem @classmethod @@ -2037,6 +2348,9 @@ def readMPS(cls, mps_file): data_model = ParseMps(mps_file) problem._from_data_model(data_model) problem.model = data_model + # Python objects match the attached file DataModel; avoid a no-op + # structure rebuild that would rewrite the model from Python. + problem._clear_stale() return problem def writeMPS(self, mps_file): @@ -2046,8 +2360,10 @@ def writeMPS(self, mps_file): -------- >>> problem.writeMPS("model.mps") """ - if self.model is None: + if self.model is None or self._stale["structure"]: self._to_data_model() + elif any(self._stale.values()): + self._refresh_data_model_values() self.model.writeMPS(mps_file) @property @@ -2108,19 +2424,17 @@ def getCSR(self): Computes and returns the CSR representation of the constraint matrix. """ - if self.constraint_csr_matrix is not None: - return self.dict_to_object(self.constraint_csr_matrix) - csr_dict = {"row_pointers": [0], "column_indices": [], "values": []} - for constr in self.constrs: - if constr.is_quadratic: - continue - csr_dict["column_indices"].extend( - list(constr.vindex_coeff_dict.keys()) - ) - csr_dict["values"].extend(list(constr.vindex_coeff_dict.values())) - csr_dict["row_pointers"].append(len(csr_dict["column_indices"])) - self.constraint_csr_matrix = csr_dict - return self.dict_to_object(csr_dict) + if self._stale["structure"] or self.constraint_csr_matrix is None: + # Rebuild typed CSR (and DataModel) so topology matches constraints. + self._to_data_model() + # Preserve the public list-valued API while keeping typed arrays + # internally for DataModel and SciPy consumers. + return self.dict_to_object( + { + key: value.tolist() if isinstance(value, np.ndarray) else value + for key, value in self.constraint_csr_matrix.items() + } + ) def getQCSR(self): """ @@ -2154,10 +2468,18 @@ def relax(self): >>> mip_problem = problem.Problem.readMPS("MIP.mps") >>> lp_problem = problem.relax() """ - self.reset_solved_values() - relaxed_problem = copy.deepcopy(self) - vars = relaxed_problem.getVariables() - for v in vars: + # DataModel wraps a Cython C++ view that cannot be deep-copied. + # Detach it for the copy; the original keeps its caches, and the + # clone rebuilds on its first solve (needed anyway after type changes). + model, self.model = self.model, None + try: + relaxed_problem = copy.deepcopy(self) + finally: + self.model = model + + # Leave the original problem's solution intact; only wipe the clone. + relaxed_problem._reset_solved_values() + for v in relaxed_problem.getVariables(): v.VariableType = CONTINUOUS return relaxed_problem @@ -2190,13 +2512,18 @@ def populate_solution(self, solution): dual_sol = None if not IsMIP: dual_sol = solution.get_dual_solution() + + # Slack is computed when Solution is built (same path as primal/dual). + slacks = solution.get_slack() + linear_row = 0 for constr in self.constrs: if constr.is_quadratic: continue if dual_sol is not None and len(dual_sol) > linear_row: constr.DualValue = dual_sol[linear_row] - constr.Slack = constr.compute_slack() + if slacks is not None: + constr.Slack = slacks[linear_row] linear_row += 1 self.solved = True @@ -2216,8 +2543,14 @@ def solve(self, settings=solver_settings.SolverSettings()): >>> problem.setObjective(x + y, sense=MAXIMIZE) >>> problem.solve() """ - if self.model is None: + if ( + self.model is None + or self.constraint_csr_matrix is None + or self._stale["structure"] + ): self._to_data_model() + else: + self._refresh_data_model_values() # Call Solver solution = solver.Solve(self.model, settings) # Post Solve diff --git a/python/cuopt/cuopt/linear_programming/solution/solution.py b/python/cuopt/cuopt/linear_programming/solution/solution.py index a1f69a6189..24a548fa19 100644 --- a/python/cuopt/cuopt/linear_programming/solution/solution.py +++ b/python/cuopt/cuopt/linear_programming/solution/solution.py @@ -74,6 +74,10 @@ class Solution: Note: Applicable to only LP The reduced cost. It contains the dual multipliers for the linear constraints. + slack : numpy.array + Slack/surplus per linear constraint in CSR row order: + ``rhs - lhs`` for ``<=``, ``lhs - rhs`` for ``>=``. Equality rows + store the residual ``rhs - lhs``. None when it was not computed. termination_status: Integer Termination status value. primal_residual: Float64 @@ -168,10 +172,12 @@ def __init__( max_variable_bound_violation=0.0, num_nodes=0, num_simplex_iterations=0, + slack=None, ): self.problem_category = problem_category self.primal_solution = primal_solution self.dual_solution = dual_solution + self.slack = slack if problem_category == ProblemCategory.LP: self.pdlp_warm_start_data = PDLPWarmStartData( current_primal_solution, @@ -262,6 +268,22 @@ def get_dual_solution(self): self.raise_if_milp_solution("get_dual_solution") return self.dual_solution + def get_slack(self): + """ + Returns the constraint slack/surplus as numpy.array with float64 type. + + For each linear constraint with ``lhs = A_i @ primal``: + + * ``<=``: ``rhs - lhs`` (slack; non-negative if feasible) + * ``>=``: ``lhs - rhs`` (surplus; non-negative if feasible) + * ``==``: ``rhs - lhs`` (residual; near zero if feasible) + + Quadratic constraints are not included. Returns None when the values + could not be computed, for example when the solver returned no primal + solution. + """ + return self.slack + def get_primal_objective(self): """ Returns the primal objective as a float64. diff --git a/python/cuopt/cuopt/linear_programming/solver/solver_wrapper.pyx b/python/cuopt/cuopt/linear_programming/solver/solver_wrapper.pyx index 1bbb101af9..7a27a87140 100644 --- a/python/cuopt/cuopt/linear_programming/solver/solver_wrapper.pyx +++ b/python/cuopt/cuopt/linear_programming/solver/solver_wrapper.pyx @@ -1,4 +1,4 @@ -# SPDX-FileCopyrightText: Copyright (c) 2023-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. # noqa +# SPDX-FileCopyrightText: Copyright (c) 2023-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. # SPDX-License-Identifier: Apache-2.0 @@ -6,6 +6,9 @@ # distutils: language = c++ # cython: embedsignature = True # cython: language_level = 3 +# cython: boundscheck=False +# cython: wraparound=False +# cython: cdivision=True from pylibraft.common.handle cimport * @@ -247,17 +250,84 @@ cdef set_solver_setting( data_model_obj.get_initial_dual_solution().shape[0] ) +cdef object _compute_slack_csr(const double[::1] rhs, + const signed char[::1] sense, + const double[::1] values, + const int[::1] indices, + const int[::1] indptr, + const double[::1] primal_solution): + """ + Classical LP slack/surplus (+ EQ residual). + + LE ('L'): rhs - lhs; GE ('G'): lhs - rhs; EQ ('E'): rhs - lhs. + + ``sense`` is signed char, not plain char: NumPy exposes an ``S1`` array as + format ``1s``, which Cython matches against signed char. Plain char is + unsigned on aarch64, so the buffer would be rejected there. + """ + cdef Py_ssize_t m = indptr.shape[0] - 1 + cdef Py_ssize_t i, k, start, end + cdef double lhs + cdef signed char s + cdef double[::1] out = np.empty(m, dtype=np.float64) + + with nogil: + for i in range(m): + start = indptr[i] + end = indptr[i + 1] + lhs = 0.0 + for k in range(start, end): + lhs = lhs + values[k] * primal_solution[indices[k]] + s = sense[i] + if s == b'G': + out[i] = lhs - rhs[i] + elif s == b'L': + out[i] = rhs[i] - lhs + else: + out[i] = rhs[i] - lhs + + return np.asarray(out) + + +cdef object _slack_from_data_model(object data_model_obj, + object primal_solution): + """Classical LE/GE slack/surplus (EQ residual) from the solved DataModel.""" + # gRPC result path has no DataModel; empty primal means no usable solution. + if data_model_obj is None or primal_solution is None: + return None + if len(primal_solution) == 0: + return None + + cdef DataModel dm = data_model_obj + offsets = dm.get_constraint_matrix_offsets() + if len(offsets) == 0: + return np.empty(0, dtype=np.float64) + + return _compute_slack_csr( + np.ascontiguousarray(dm.get_constraint_bounds(), dtype=np.float64), + np.ascontiguousarray(dm.get_row_types(), dtype="S1"), + np.ascontiguousarray(dm.get_constraint_matrix_values(), dtype=np.float64), + np.ascontiguousarray(dm.get_constraint_matrix_indices(), dtype=np.int32), + np.ascontiguousarray(offsets, dtype=np.int32), + np.ascontiguousarray(primal_solution, dtype=np.float64), + ) + + cdef create_solution(unique_ptr[solver_ret_t] sol_ret_ptr, DataModel data_model_obj, is_batch=False): return create_solution_with_names( - move(sol_ret_ptr), data_model_obj.get_variable_names(), is_batch + move(sol_ret_ptr), + data_model_obj.get_variable_names(), + is_batch, + data_model_obj, ) cdef create_solution_with_names(unique_ptr[solver_ret_t] sol_ret_ptr, object variable_names, - bint is_batch=False): + bint is_batch=False, + object data_model_obj=None): from cuopt.linear_programming.solution.solution import Solution @@ -294,7 +364,8 @@ cdef create_solution_with_names(unique_ptr[solver_ret_t] sol_ret_ptr, max_int_violation=mip_ptr.max_int_violation_, max_constraint_violation=mip_ptr.max_constraint_violation_, num_nodes=mip_ptr.nodes_, - num_simplex_iterations=mip_ptr.simplex_iterations_ + num_simplex_iterations=mip_ptr.simplex_iterations_, + slack=_slack_from_data_model(data_model_obj, solution), ) else: @@ -417,6 +488,7 @@ cdef create_solution_with_names(unique_ptr[solver_ret_t] sol_ret_ptr, lp_ptr.gap_, lp_ptr.nb_iterations_, lp_ptr.solved_by_, + slack=_slack_from_data_model(data_model_obj, primal_solution), ) else: return Solution( @@ -436,6 +508,7 @@ cdef create_solution_with_names(unique_ptr[solver_ret_t] sol_ret_ptr, gap=lp_ptr.gap_, nb_iterations=lp_ptr.nb_iterations_, solved_by=lp_ptr.solved_by_, + slack=_slack_from_data_model(data_model_obj, primal_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 964809553f..1984f9ab72 100644 --- a/python/cuopt/cuopt/tests/linear_programming/test_python_API.py +++ b/python/cuopt/cuopt/tests/linear_programming/test_python_API.py @@ -109,7 +109,7 @@ def test_model(): assert csr.column_indices == expected_column_indices assert csr.values == expected_values - expected_slack = [-6, 0] + expected_slack = [6, 0] expected_var_values = [36, 41] for i, var in enumerate(prob.getVariables()): @@ -153,11 +153,18 @@ def test_model(): def test_constraint_duplicate_terms_slack(): """Merged coeffs in vindex_coeff_dict must not be double-counted in slack.""" prob = Problem() - x = prob.addVariable() + x = prob.addVariable(lb=0.0, ub=10.0, obj=0.0) c = prob.addConstraint(5 * x + 7 * x <= 18) assert c.getCoefficient(x) == 12 - x.Value = 1.0 - assert c.compute_slack() == pytest.approx(6.0) + # Feasible point with slack 6 under classical LE slack: rhs - A@x. + settings = SolverSettings() + settings.set_parameter("time_limit", 10) + # Fix x=1 via bounds so populate_solution computes a known slack. + x.LB = 1.0 + x.UB = 1.0 + prob.setObjective(0 * x, sense=MINIMIZE) + prob.solve(settings) + assert c.Slack == pytest.approx(6.0) def test_semi_continuous_variable(): diff --git a/python/cuopt/cuopt/tests/socp/test_socp.py b/python/cuopt/cuopt/tests/socp/test_socp.py index 8235645476..1f44b69d3c 100644 --- a/python/cuopt/cuopt/tests/socp/test_socp.py +++ b/python/cuopt/cuopt/tests/socp/test_socp.py @@ -14,6 +14,7 @@ import numpy as np import pytest +from cuopt.linear_programming import Read from cuopt.linear_programming.problem import EQ, GE, LE, MAXIMIZE, Problem from cuopt.linear_programming.solver.solver_parameters import CUOPT_METHOD from cuopt.linear_programming.solver_settings import ( @@ -119,11 +120,10 @@ def _assert_feasible(problem: Problem) -> None: _quadratic_constraint_violation(constr, variables) <= FEAS_TOL ) continue - slack = constr.compute_slack() - if constr.Sense == LE: + # Classical slack/surplus from populate_solution (non-negative if feasible). + slack = constr.Slack + if constr.Sense in (LE, GE): assert slack >= -FEAS_TOL - elif constr.Sense == GE: - assert slack <= FEAS_TOL else: assert constr.Sense == EQ assert slack == pytest.approx(0.0, abs=FEAS_TOL) @@ -214,3 +214,88 @@ def test_maximize_with_quadratic_constraint(): assert prob_max.ObjValue == pytest.approx(2.0, abs=OBJ_TOL) assert x.Value == pytest.approx(1.0, abs=PRIMAL_TOL) assert y.Value == pytest.approx(1.0, abs=PRIMAL_TOL) + + +# Same model as test_maximize_with_quadratic_constraint written as MPS (as a +# minimization). QC0 is the binding row: without it the optimum is -10. +QC_MPS = """NAME QCREAD +ROWS + N OBJ + L LIN0 + L QC0 +COLUMNS + x OBJ -1 + x LIN0 1 + y OBJ -1 + y LIN0 1 +RHS + RHS1 LIN0 10 + RHS1 QC0 6 +QCMATRIX QC0 + x x 2 + x y 1 + y x 1 + y y 2 +BOUNDS + MI BND x + MI BND y +ENDATA +""" + + +def _write_qc_mps(tmp_path) -> str: + path = tmp_path / "qc_read.mps" + path.write_text(QC_MPS) + return str(path) + + +def test_read_keeps_quadratic_constraints(tmp_path): + """QCMATRIX rows parsed by read() must reach the solver.""" + path = _write_qc_mps(tmp_path) + problem = Problem.read(path) + + assert problem.NumConstraints == 2 + quad = problem.getQuadraticConstraints() + assert len(quad) == 1 + + # Each row must mirror the bundle the reader parsed, in whatever form the + # reader normalized it to. + bundle = Read(path).get_quadratic_constraints()[0] + assert quad[0].ConstraintName == bundle["constraint_row_name"] + assert quad[0].Sense == bundle["constraint_row_type"] + assert quad[0].rhs_value == pytest.approx(bundle["rhs_value"]) + for field in ("rows", "cols", "vals", "linear_indices", "linear_values"): + np.testing.assert_allclose(getattr(quad[0], field), bundle[field]) + + # The file DataModel is the model to solve; rebuilding it from Python is + # what used to drop the quadratic rows. + model = problem.model + solution = _solve(problem) + assert problem.model is model + + _assert_solution_on_original_model(problem, solution) + _assert_feasible(problem) + + x, y = problem.getVariables() + assert problem.ObjValue == pytest.approx(-2.0, abs=OBJ_TOL) + assert x.Value == pytest.approx(1.0, abs=PRIMAL_TOL) + assert y.Value == pytest.approx(1.0, abs=PRIMAL_TOL) + + +def test_read_then_rebuild_keeps_quadratic_constraints(tmp_path): + """A rebuild after read() must re-emit the quadratic rows.""" + problem = Problem.read(_write_qc_mps(tmp_path)) + # update() drops the cached CSR, so writeMPS rebuilds the DataModel from + # the Python objects instead of reusing the one that was read. + problem.getConstraint(0).RHS = 9.0 + problem.update() + + round_trip = tmp_path / "round_trip.mps" + problem.writeMPS(str(round_trip)) + assert "QCMATRIX" in round_trip.read_text() + + reread = Problem.read(str(round_trip)) + assert len(reread.getQuadraticConstraints()) == 1 + _solve(reread) + _assert_feasible(reread) + assert reread.ObjValue == pytest.approx(-2.0, abs=OBJ_TOL)