diff --git a/.gitignore b/.gitignore index 6794dcf..2a80fd4 100644 --- a/.gitignore +++ b/.gitignore @@ -58,12 +58,18 @@ dkms.conf # ignored files build/* test/* +!test/*.py /.vscode /.venv /_b *.whl *.pyc +# coverage / pytest artifacts +.coverage +coverage.xml +htmlcov/ +.pytest_cache/ *.txt *.sh \ No newline at end of file diff --git a/pyproject.toml b/pyproject.toml index 8ac6244..62b0535 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -52,7 +52,7 @@ test = [ ] [tool.pytest.ini_options] -testpaths = ["tests"] +testpaths = ["test"] addopts = """ -q -ra diff --git a/pytest.ini b/pytest.ini deleted file mode 100644 index d2a5db1..0000000 --- a/pytest.ini +++ /dev/null @@ -1,2 +0,0 @@ -[pytest] -addopts = -q \ No newline at end of file diff --git a/python/README.md b/python/README.md index 898349a..1c608e9 100644 --- a/python/README.md +++ b/python/README.md @@ -75,7 +75,7 @@ m.setParams(OutputFlag=True, OptimalityTol=1e-8) m.optimize() # Retrieve results -print("Status:", m.Status) +print("Status:", m.StatusName) print("Objective:", m.ObjVal) print("Primal solution:", m.X) print("Dual solution:", m.Pi) @@ -98,7 +98,7 @@ $$ - **constraint_matrix** (`A`): Coefficient matrix for the constraints. Both dense (`numpy.ndarray`) and sparse (`scipy.sparse.csr_matrix`) inputs are supported. Internally stored in double precision (`float64`). - **constraint_lower_bound** (`l`): Lower bounds for each constraint. Use `-np.inf` or `None` for no lower bound. - **constraint_upper_bound** (`u`): Upper bounds for each constraint. Use `+np.inf` or `None` for no upper bound. -- **variable_lower_bound** (`lb`, optional): Lower bounds for the decision variables. Defaults to `0` for all variables if not provided. +- **variable_lower_bound** (`lb`, optional): Lower bounds for the decision variables. Defaults to `-np.inf` for all variables if not provided. - **variable_upper_bound** (`ub`, optional): Upper bounds for the decision variables. Defaults to `+np.inf` for all variables if not provided. - **objective_constant** (`c0`, optional): Constant offset in the objective function. Defaults to `0.0`. @@ -125,7 +125,7 @@ m = Model(objective_vector=c, ### Reading from MPS Files -A `Model` can also be created directly from an MPS file (plain or gzip-compressed) with `cupdlpx.read`, similar to `gurobipy.read`: +A `Model` can also be created directly from an MPS file (plain or gzip-compressed) with `cupdlpx.read`: ```python import cupdlpx @@ -174,7 +174,7 @@ Below is a list of commonly used parameters, their internal keys, and descriptio | `ReflectionCoeff` | `reflection_coefficient` | float | `1.0` | Reflection coefficient. | | `SVMaxIter` | `sv_max_iter` | int | 5000 | Maximum number of iterations for the power method | | `SVTol`| `sv_tol` | float | `1e-4` | Termination tolerance for the power method | -| `Presolve`| `presolve` | float | `True` | Whether to use presolve. | +| `Presolve`| `presolve` | bool | `True` | Whether to use presolve. | | `FeasibilityPolishing` | `feasibility_polishing` | bool | `False` | Run feasibility polishing process.| | `FeasibilityPolishingTol` | `eps_feas_polish_relative` | float | `1e-6` | Relative tolerance for primal/dual residual. | @@ -191,18 +191,21 @@ m.setParams(TimeLimit=300, FeasibilityTol=1e-6) # Method 3: attribute-style access m.Params.TimeLimit = 300 m.Params.FeasibilityTol = 1e-6 + +# Reset all parameters to backend defaults +m.resetParams() ``` ## Solution Attributes -After calling `m.optimize()`, the solver stores results in a set of read-only attributes. These attributes provide access to primal/dual solutions, objective values, residuals, and runtime statistics. +After calling `m.optimize()`, the solver stores results in a set of read-only attributes. `optimize()` returns the model itself, so chained access like `m.optimize().Status` is also supported. These attributes provide access to primal/dual solutions, objective values, residuals, and runtime statistics. ### Attribute Reference | Attribute | Type | Description | |---|---|---| -| `Status` | str | Human-readable solver status (`"OPTIMAL"`, `"INFEASIBLE"`, `"UNBOUNDED"`, `"TIME_LIMIT"`, etc.). | -| `StatusCode` | int | Numeric status code (`OPTIMAL=1`, `INFEASIBLE=2`, `UNBOUNDED=3`, `ITERATION_LIMIT=4`, `TIME_LIMIT=5`, `UNSPECIFIED=-1`). | +| `Status` | int | Integer termination status code; compare against `cupdlpx.PDLP` constants: `OPTIMAL=0`, `PRIMAL_INFEASIBLE=1`, `DUAL_INFEASIBLE=2`, `TIME_LIMIT=3`, `ITERATION_LIMIT=4`, `INFEASIBLE_OR_UNBOUNDED=5`, `FEAS_POLISH_SUCCESS=6`, `UNSPECIFIED=-1`. | +| `StatusName` | str | Human-readable status name, e.g. `"OPTIMAL"`, `"PRIMAL_INFEASIBLE"`. | | `ObjVal` | float | Primal objective value at termination (sign-adjusted according to `ModelSense`). | | `DualObj` | float | Dual objective value at termination. | | `Gap` | float | Absolute primal-dual gap. | @@ -225,7 +228,7 @@ All solution-related information can then be queried directly from the `Model` o ```python m.optimize() -print("Status:", m.Status, "(code:", m.StatusCode, ")") +print("Status:", m.StatusName, "(code:", m.Status, ")") print("Primal objective:", m.ObjVal) print("Dual objective:", m.DualObj) print("Relative gap:", m.RelGap) @@ -268,7 +271,7 @@ m.setWarmStart(primal=x_init) m.setWarmStart(dual=pi_init) ``` -If the warm-start vectors have incorrect dimensions, the solver automatically falls back to a cold start and issues a warning. +If the warm-start vectors have incorrect dimensions, `setWarmStart` raises a `ValueError`. Omitting an argument leaves that side unchanged; passing `None` clears it. To clear existing warm-start values: diff --git a/python/cupdlpx/PDLP.py b/python/cupdlpx/PDLP.py index 0dfa1f0..e990adb 100644 --- a/python/cupdlpx/PDLP.py +++ b/python/cupdlpx/PDLP.py @@ -12,17 +12,21 @@ # See the License for the specific language governing permissions and # limitations under the License. -# Objective +"""Objective-sense and termination-status constants, plus parameter-name aliases.""" + +# Objective MINIMIZE = 1 MAXIMIZE = -1 -# Status codes -OPTIMAL = 0 -PRIMAL_INFEASIBLE = 1 -DUAL_INFEASIBLE = 2 -TIME_LIMIT = 3 -ITERATION_LIMIT = 4 -UNSPECIFIED = -1 +# Status codes (must match status_to_code in python_bindings/_core_bindings.cpp) +OPTIMAL = 0 +PRIMAL_INFEASIBLE = 1 +DUAL_INFEASIBLE = 2 +TIME_LIMIT = 3 +ITERATION_LIMIT = 4 +INFEASIBLE_OR_UNBOUNDED = 5 +FEAS_POLISH_SUCCESS = 6 +UNSPECIFIED = -1 # parameter name alias diff --git a/python/cupdlpx/__init__.py b/python/cupdlpx/__init__.py index 00ea958..fed1f83 100644 --- a/python/cupdlpx/__init__.py +++ b/python/cupdlpx/__init__.py @@ -12,11 +12,13 @@ # See the License for the specific language governing permissions and # limitations under the License. +"""cuPDLPx: Python bindings for the GPU-accelerated first-order LP solver.""" + import os import platform # Windows only: register CUDA bin for dependent DLL loading. -if platform.system() == "Windows": +if platform.system() == "Windows": # pragma: no cover cuda_path = os.environ.get("CUDA_PATH") if cuda_path: bin_path = os.path.join(cuda_path, "bin") @@ -26,12 +28,12 @@ from .model import Model, read from . import PDLP -__all__ = ["Model", "read"] - # versioning from importlib.metadata import version, PackageNotFoundError # get version from package metadata (toml file) try: __version__ = version("cupdlpx") -except PackageNotFoundError: +except PackageNotFoundError: # pragma: no cover __version__ = "0.0.0" + +__all__ = ["Model", "PDLP", "read", "__version__"] diff --git a/python/cupdlpx/_core.py b/python/cupdlpx/_core.py index dd87cfb..cdd7c95 100644 --- a/python/cupdlpx/_core.py +++ b/python/cupdlpx/_core.py @@ -12,4 +12,6 @@ # See the License for the specific language governing permissions and # limitations under the License. +"""Thin re-export of the compiled cuPDLPx core extension (_cupdlpx_core).""" + from ._cupdlpx_core import solve_once, get_default_params, read_mps \ No newline at end of file diff --git a/python/cupdlpx/model.py b/python/cupdlpx/model.py index 87239f6..82d2a7e 100644 --- a/python/cupdlpx/model.py +++ b/python/cupdlpx/model.py @@ -12,9 +12,10 @@ # See the License for the specific language governing permissions and # limitations under the License. +"""Model interface for the cuPDLPx LP solver.""" + from __future__ import annotations import os -import warnings from typing import Any, Optional, Union import numpy as np @@ -26,47 +27,144 @@ # array-like type ArrayLike = Union[np.ndarray, list, tuple] +# sentinel for "argument not provided" (distinct from None, which means "clear") +_UNSET = object() + +_BOOL_PARAMS = frozenset( + { + "verbose", + "has_pock_chambolle_alpha", + "bound_objective_rescaling", + "feasibility_polishing", + "presolve", + } +) +_INT_PARAMS = frozenset( + { + "termination_evaluation_frequency", + "iteration_limit", + "l_inf_ruiz_iterations", + "sv_max_iter", + } +) +_POSITIVE_INT_PARAMS = frozenset({"termination_evaluation_frequency", "sv_max_iter"}) +_FLOAT_PARAMS = frozenset( + { + "eps_optimal_relative", + "eps_feasible_relative", + "time_sec_limit", + "pock_chambolle_alpha", + "artificial_restart_threshold", + "sufficient_reduction_for_restart", + "necessary_reduction_for_restart", + "k_p", + "reflection_coefficient", + "eps_feas_polish_relative", + "sv_tol", + "matrix_zero_tol", + } +) +_POSITIVE_FLOAT_PARAMS = frozenset( + { + "eps_optimal_relative", + "eps_feasible_relative", + "eps_feas_polish_relative", + "sv_tol", + } +) +_NONNEGATIVE_FLOAT_PARAMS = frozenset({"time_sec_limit", "matrix_zero_tol"}) +_STRING_PARAMS = frozenset({"optimality_norm"}) + +# int params are stored as C int32 on the backend +_INT32_MAX = np.iinfo(np.int32).max + +# every backend param must fall in one typed set, else it goes unvalidated +_CLASSIFIED_PARAMS = _BOOL_PARAMS | _INT_PARAMS | _FLOAT_PARAMS | _STRING_PARAMS + +# refinement sets must be subsets of their base type sets +assert _POSITIVE_INT_PARAMS <= _INT_PARAMS +assert _POSITIVE_FLOAT_PARAMS <= _FLOAT_PARAMS +assert _NONNEGATIVE_FLOAT_PARAMS <= _FLOAT_PARAMS + def _as_dense_f64_c(a: ArrayLike) -> np.ndarray: """ - Convert input to a C-contiguous numpy array of float64. + Convert input to an owned C-contiguous numpy array of float64. """ - arr = np.asarray(a, dtype=np.float64) - # ensure C-contiguous - if not arr.flags.c_contiguous: - arr = np.ascontiguousarray(arr, dtype=np.float64) + return np.array(a, dtype=np.float64, order="C", copy=True) + +def _readonly_array(arr: np.ndarray) -> np.ndarray: + arr.setflags(write=False) return arr -def _as_csr_f64_i32(A: sp.spmatrix) -> sp.csr_matrix: +def _readonly_matrix(A): + if sp.issparse(A): + A.data.setflags(write=False) + A.indices.setflags(write=False) + A.indptr.setflags(write=False) + else: + A.setflags(write=False) + return A + +def _require_finite(name: str, arr: np.ndarray) -> None: + if not np.all(np.isfinite(arr)): + raise ValueError(f"{name} must contain only finite values") + +def _require_no_nan(name: str, arr: np.ndarray) -> None: + if np.any(np.isnan(arr)): + raise ValueError(f"{name} must not contain NaN") + +def _check_bounds(lower: Optional[np.ndarray], upper: Optional[np.ndarray], name: str) -> None: + if lower is not None and upper is not None and np.any(lower > upper): + raise ValueError(f"{name}: lower bounds must be <= upper bounds") + +def _as_csr_f64_i32(A) -> sp.csr_matrix: """ - Convert input sparse matrix to CSR format with float64 values and int32 indices. + Convert input sparse matrix/array to CSR format with float64 values and + int32 indices. Never mutates the caller's matrix. """ - csr = A.tocsr().astype(np.float64, copy=False) - # force int32 indices (common C/CUDA req) + csr = A.tocsr() + if csr is A: + csr = csr.copy() + if csr.dtype != np.float64: + csr = csr.astype(np.float64) + _require_finite("constraint_matrix", csr.data) + # merge duplicate (row, col) entries and sort indices within each row + csr.sum_duplicates() + csr.sort_indices() + # force int32 indices (common C/CUDA req); int32 arrays can't overflow if csr.indptr.dtype != np.int32: - csr.indptr = csr.indptr.astype(np.int32, copy=True) + # indptr is non-decreasing, so its last entry (== nnz) is the maximum + if csr.indptr[-1] > _INT32_MAX: + raise OverflowError("constraint_matrix CSR indptr exceeds int32 range") + csr.indptr = csr.indptr.astype(np.int32, copy=False) if csr.indices.dtype != np.int32: - csr.indices = csr.indices.astype(np.int32, copy=True) - csr.sort_indices() - return csr + if csr.indices.size and csr.indices.max() > _INT32_MAX: + raise OverflowError("constraint_matrix CSR indices exceed int32 range") + csr.indices = csr.indices.astype(np.int32, copy=False) + return _readonly_matrix(csr) def read(filename: Union[str, os.PathLike]) -> "Model": """ Read a linear program from an MPS file (plain or gzip-compressed) and - return a Model, similar to gurobipy.read. + return a Model. Parameters: - filename: Path to a .mps or .mps.gz file. """ + # normalize path and check existence filename = os.fspath(filename) if not os.path.isfile(filename): raise FileNotFoundError(f"No such MPS file: {filename}") + # parse the MPS file into raw problem data data = read_mps(str(filename)) + # rebuild the constraint matrix in CSR form m = int(data["num_constraints"]) n = int(data["num_variables"]) A = sp.csr_matrix( (data["values"], data["col_ind"], data["row_ptr"]), shape=(m, n) ) + # assemble the model model = Model( objective_vector=data["objective_vector"], constraint_matrix=A, @@ -76,6 +174,7 @@ def read(filename: Union[str, os.PathLike]) -> "Model": variable_upper_bound=data["variable_upper_bound"], objective_constant=data["objective_constant"], ) + # set objective sense model.ModelSense = PDLP.MAXIMIZE if data["maximize"] else PDLP.MINIMIZE return model @@ -105,6 +204,25 @@ def __setitem__(self, name: str, value): def keys(self): return self._m._params.keys() + def values(self): + return self._m._params.values() + + def items(self): + return self._m._params.items() + + def __contains__(self, name: str): + key = PDLP._PARAM_ALIAS.get(name, name) + return key in self._m._params + + def __iter__(self): + return iter(self._m._params) + + def __len__(self): + return len(self._m._params) + + def __repr__(self): + return f"ParamsView({dict(self._m._params)!r})" + class Model: """ @@ -125,14 +243,16 @@ def __init__( Parameters: - objective_vector: Coefficients of the objective function. - - constraint_matrix: Coefficients of the constraints. - - lower_bounds: Lower bounds for the decision variables. - - upper_bounds: Upper bounds for the decision variables. - - constraint_lower_bounds: Lower bounds for the constraints. - - constraint_upper_bounds: Upper bounds for the constraints. + - constraint_matrix: Constraint coefficient matrix (2D dense or scipy.sparse). + - constraint_lower_bound: Lower bounds for the constraints. + - constraint_upper_bound: Upper bounds for the constraints. + - variable_lower_bound: Lower bounds for the decision variables (default -inf). + - variable_upper_bound: Upper bounds for the decision variables (default +inf). - objective_constant: Constant term in the objective function. - - model_sense: PDLP.MINIMIZE or PDLP.MAXIMIZE. - If variable bounds are not provided, they default to -inf and +inf respectively. + + The objective sense defaults to PDLP.MINIMIZE; set model.ModelSense to + PDLP.MAXIMIZE to maximize. Constraint bounds may be None, meaning -inf + (lower) or +inf (upper). """ # problem dimensions if not hasattr(constraint_matrix, "shape") or len(constraint_matrix.shape) != 2: @@ -140,23 +260,30 @@ def __init__( m, n = constraint_matrix.shape self.num_vars = int(n) self.num_constrs = int(m) - # sense - self.ModelSense = PDLP.MINIMIZE + # model data storage (populated by the setters below; exposed via properties) + self._A = None + self._c: Optional[np.ndarray] = None + self._c0: float = 0.0 + self._lb: Optional[np.ndarray] = None + self._ub: Optional[np.ndarray] = None + self._constr_lb: Optional[np.ndarray] = None + self._constr_ub: Optional[np.ndarray] = None + # objective sense (default minimize) + self._model_sense = PDLP.MINIMIZE # always start from backend defaults PDLP params - self._params: dict[str, Any] = dict(get_default_params()) + self._default_params: dict[str, Any] = dict(get_default_params()) + self._params: dict[str, Any] = dict(self._default_params) + # canonical set of backend parameter keys, used to reject typos in setParam + self._valid_param_keys = frozenset(self._params) + # fail loudly if the backend exposes a param the typed sets don't cover + unclassified = self._valid_param_keys - _CLASSIFIED_PARAMS + if unclassified: + raise RuntimeError(f"Unclassified solver parameters (update model.py): {sorted(unclassified)}") self.Params = _ParamsView(self) - # set coefficients and bounds - self.setObjectiveVector(objective_vector) - self.setObjectiveConstant(objective_constant) - self.setConstraintMatrix(constraint_matrix) - self.setConstraintLowerBound(constraint_lower_bound) - self.setConstraintUpperBound(constraint_upper_bound) - self.setVariableLowerBound(variable_lower_bound) - self.setVariableUpperBound(variable_upper_bound) # initialize warm start values self._primal_start: Optional[np.ndarray] = None # warm start primal solution self._dual_start: Optional[np.ndarray] = None # warm start dual solution - # initialize solution attributes + # initialize solution attributes before the setters below self._x: Optional[np.ndarray] = None # primal solution self._y: Optional[np.ndarray] = None # dual solution self._rc: Optional[np.ndarray] = None # reduced costs @@ -164,8 +291,8 @@ def __init__( self._dualobj: Optional[float] = None # dual objective value self._gap: Optional[float] = None # primal-dual gap self._rel_gap: Optional[float] = None # relative gap - self._status: Optional[str] = None # solution status - self._status_code: Optional[int] = None # solution status code + self._status_name: Optional[str] = None # solution status name (str) + self._status_code: Optional[int] = None # solution status code (int) self._iter: Optional[int] = None # number of iterations self._runtime: Optional[float] = None # runtime self._rescale_time: Optional[float] = None # rescale time @@ -175,18 +302,27 @@ def __init__( self._max_d_ray: Optional[float] = None # maximum dual ray self._p_ray_lin_obj: Optional[float] = None # primal ray linear objective self._d_ray_obj: Optional[float] = None # dual ray objective + # set coefficients and bounds + self.setObjectiveVector(objective_vector) + self.setObjectiveConstant(objective_constant) + self.setConstraintMatrix(constraint_matrix) + self.setConstraintLowerBound(constraint_lower_bound) + self.setConstraintUpperBound(constraint_upper_bound) + self.setVariableLowerBound(variable_lower_bound) + self.setVariableUpperBound(variable_upper_bound) + self._validate_bounds() def setObjectiveVector(self, c: ArrayLike) -> None: """ Overwrite objective vector c. """ - # store as float64 - self.c = _as_dense_f64_c(c) - # check dimensions - if self.c.ndim != 1: - raise ValueError(f"setObjectiveVector: c must be 1D, got shape {self.c.shape}") - if self.c.size != self.num_vars: - raise ValueError(f"setObjectiveVector: length {self.c.size} != self.num_vars ({self.num_vars})") + c_arr = _as_dense_f64_c(c) + if c_arr.ndim != 1: + raise ValueError(f"setObjectiveVector: c must be 1D, got shape {c_arr.shape}") + if c_arr.size != self.num_vars: + raise ValueError(f"setObjectiveVector: length {c_arr.size} != self.num_vars ({self.num_vars})") + _require_finite("objective_vector", c_arr) + self._c = _readonly_array(c_arr) # clear cached solution self._clear_solution_cache() @@ -195,47 +331,53 @@ def setObjectiveConstant(self, c0: float) -> None: Overwrite objective constant term. Minimal check: convert to float. """ - self.c0 = float(c0) + c0 = float(c0) + if not np.isfinite(c0): + raise ValueError("objective_constant must be finite") + self._c0 = c0 # clear cached solution self._clear_solution_cache() - def setConstraintMatrix(self, A_like: ArrayLike) -> None: + def setConstraintMatrix(self, A_like: Union[np.ndarray, sp.spmatrix]) -> None: """ Overwrite constraint matrix A. """ - if not isinstance(A_like, (np.ndarray, sp.spmatrix)): - raise TypeError("setConstraintMatrix: A must be a numpy.ndarray or scipy.sparse matrix") + if not (sp.issparse(A_like) or isinstance(A_like, np.ndarray)): + raise TypeError("setConstraintMatrix: A must be a numpy.ndarray or scipy.sparse matrix/array") if len(A_like.shape) != 2: raise ValueError(f"setConstraintMatrix: A must be 2D, got shape {A_like.shape}") if A_like.shape[1] != self.num_vars: raise ValueError(f"setConstraintMatrix: A shape {A_like.shape} does not match number of variables ({self.num_vars})") - # store as float64 + # convert to backend layout (do not mutate self until all checks pass) if sp.issparse(A_like): - self.A = _as_csr_f64_i32(A_like) + A = _as_csr_f64_i32(A_like) else: - self.A = _as_dense_f64_c(A_like) - # problem dimensions - if not hasattr(self.A, "shape") or len(self.A.shape) != 2: - raise ValueError("constraint_matrix must be a 2D numpy.ndarray or scipy.sparse matrix.") - m, _ = self.A.shape - self.num_constrs = int(m) - # check constraint bounds - l = getattr(self, "constr_lb", None) + A = _as_dense_f64_c(A_like) + _require_finite("constraint_matrix", A) + m = int(A.shape[0]) + # validate existing constraint bounds against the new row count before committing + l = self._constr_lb if l is not None: - l = np.asarray(l, dtype=np.float64).ravel() - if l.size != self.num_constrs: + n_l = np.asarray(l).ravel().size + if n_l != m: raise ValueError( - f"setConstraintMatrix: constraint_lower_bound length {l.size} != rows {self.num_constrs}. " + f"setConstraintMatrix: constraint_lower_bound length {n_l} != rows {m}. " f"Call setConstraintLowerBound(...) to update it." ) - u = getattr(self, "constr_ub", None) + u = self._constr_ub if u is not None: - u = np.asarray(u, dtype=np.float64).ravel() - if u.size != self.num_constrs: + n_u = np.asarray(u).ravel().size + if n_u != m: raise ValueError( - f"setConstraintMatrix: constraint_upper_bound length {u.size} != rows {self.num_constrs}. " + f"setConstraintMatrix: constraint_upper_bound length {n_u} != rows {m}. " f"Call setConstraintUpperBound(...) to update it." - ) + ) + # commit + self._A = _readonly_matrix(A) + self.num_constrs = m + # drop a dual warm start that no longer matches the row count + if self._dual_start is not None and self._dual_start.size != m: + self._dual_start = None # clear cached solution self._clear_solution_cache() @@ -245,17 +387,17 @@ def setConstraintLowerBound(self, constr_lb: Optional[ArrayLike]) -> None: """ # check if the input is None if constr_lb is None: - self.constr_lb = None + self._constr_lb = None # clear cached solution self._clear_solution_cache() return - # convert to numpy array - constr_lb = _as_dense_f64_c(constr_lb).ravel() - if constr_lb.size != self.num_constrs: + constr_lb_arr = _as_dense_f64_c(constr_lb).ravel() + if constr_lb_arr.size != self.num_constrs: raise ValueError( - f"setConstraintLowerBound: length {constr_lb.size} != self.num_constrs ({self.num_constrs})" + f"setConstraintLowerBound: length {constr_lb_arr.size} != self.num_constrs ({self.num_constrs})" ) - self.constr_lb = constr_lb + _require_no_nan("constraint_lower_bound", constr_lb_arr) + self._constr_lb = _readonly_array(constr_lb_arr) # clear cached solution self._clear_solution_cache() @@ -265,17 +407,17 @@ def setConstraintUpperBound(self, constr_ub: Optional[ArrayLike]) -> None: """ # check if the input is None if constr_ub is None: - self.constr_ub = None + self._constr_ub = None # clear cached solution self._clear_solution_cache() return - # convert to numpy array - constr_ub = _as_dense_f64_c(constr_ub).ravel() - if constr_ub.size != self.num_constrs: + constr_ub_arr = _as_dense_f64_c(constr_ub).ravel() + if constr_ub_arr.size != self.num_constrs: raise ValueError( - f"setConstraintUpperBound: length {constr_ub.size} != self.num_constrs ({self.num_constrs})" + f"setConstraintUpperBound: length {constr_ub_arr.size} != self.num_constrs ({self.num_constrs})" ) - self.constr_ub = constr_ub + _require_no_nan("constraint_upper_bound", constr_ub_arr) + self._constr_ub = _readonly_array(constr_ub_arr) # clear cached solution self._clear_solution_cache() @@ -285,17 +427,17 @@ def setVariableLowerBound(self, lb: Optional[ArrayLike]) -> None: """ # check if the input is None if lb is None: - self.lb = None + self._lb = None # clear cached solution self._clear_solution_cache() return - # convert to numpy array - lb = _as_dense_f64_c(lb).ravel() - if lb.size != self.num_vars: + lb_arr = _as_dense_f64_c(lb).ravel() + if lb_arr.size != self.num_vars: raise ValueError( - f"setVariableLowerBound: length {lb.size} != self.num_vars ({self.num_vars})" + f"setVariableLowerBound: length {lb_arr.size} != self.num_vars ({self.num_vars})" ) - self.lb = lb + _require_no_nan("variable_lower_bound", lb_arr) + self._lb = _readonly_array(lb_arr) # clear cached solution self._clear_solution_cache() @@ -305,50 +447,56 @@ def setVariableUpperBound(self, ub: Optional[ArrayLike]) -> None: """ # check if the input is None if ub is None: - self.ub = None + self._ub = None # clear cached solution self._clear_solution_cache() return - # convert to numpy array - ub = _as_dense_f64_c(ub).ravel() - if ub.size != self.num_vars: + ub_arr = _as_dense_f64_c(ub).ravel() + if ub_arr.size != self.num_vars: raise ValueError( - f"setVariableUpperBound: length {ub.size} != self.num_vars ({self.num_vars})" + f"setVariableUpperBound: length {ub_arr.size} != self.num_vars ({self.num_vars})" ) - self.ub = ub + _require_no_nan("variable_upper_bound", ub_arr) + self._ub = _readonly_array(ub_arr) # clear cached solution self._clear_solution_cache() - def setWarmStart(self, primal: Optional[ArrayLike] = None, dual: Optional[ArrayLike] = None) -> None: + def setWarmStart(self, primal: Optional[ArrayLike] = _UNSET, dual: Optional[ArrayLike] = _UNSET) -> None: """ - Set warm start values for primal and/or dual solutions. + Set warm start values for the primal and/or dual solutions. + + For each of primal/dual: pass an array to set it, None to clear it, or + omit the argument to leave the current value unchanged. Raises + ValueError on a size mismatch. """ - # set primal warm start - if primal is not None: - primal_arr = _as_dense_f64_c(primal).ravel() - if primal_arr.size == self.num_vars: # otherwise default to None - self._primal_start = primal_arr + next_primal = self._primal_start + next_dual = self._dual_start + # primal warm start + if primal is not _UNSET: + if primal is None: + next_primal = None else: - warnings.warn( - f"Warm start primal size mismatch (expected {self.num_vars}, got {primal_arr.size}).", - RuntimeWarning - ) - # clear primal warm start - else: - self._primal_start = None - # set dual warm start - if dual is not None: - dual_arr = _as_dense_f64_c(dual).ravel() - if dual_arr.size == self.num_constrs: # otherwise default to None - self._dual_start = dual_arr + primal_arr = _as_dense_f64_c(primal).ravel() + if primal_arr.size != self.num_vars: + raise ValueError( + f"setWarmStart: primal size mismatch (expected {self.num_vars}, got {primal_arr.size})." + ) + _require_finite("primal warm start", primal_arr) + next_primal = _readonly_array(primal_arr) + # dual warm start + if dual is not _UNSET: + if dual is None: + next_dual = None else: - warnings.warn( - f"Warm start dual size mismatch (expected {self.num_constrs}, got {dual_arr.size}).", - RuntimeWarning - ) - # clear dual warm start - else: - self._dual_start = None + dual_arr = _as_dense_f64_c(dual).ravel() + if dual_arr.size != self.num_constrs: + raise ValueError( + f"setWarmStart: dual size mismatch (expected {self.num_constrs}, got {dual_arr.size})." + ) + _require_finite("dual warm start", dual_arr) + next_dual = _readonly_array(dual_arr) + self._primal_start = next_primal + self._dual_start = next_dual def clearWarmStart(self) -> None: """ @@ -356,26 +504,104 @@ def clearWarmStart(self) -> None: """ self.setWarmStart(primal=None, dual=None) + def _resolve_param_key(self, name: str) -> str: + """ + Map a user-facing parameter name (alias or backend key) to its backend + key, raising KeyError for unknown names instead of silently accepting them. + """ + # map alias to backend key + key = PDLP._PARAM_ALIAS.get(name, name) + # reject unknown names + if key not in self._valid_param_keys: + valid = sorted(PDLP._PARAM_ALIAS.keys()) + sorted(self._valid_param_keys) + raise KeyError(f"Unknown parameter '{name}'. Valid names: {valid}") + return key + + def _validate_param_value(self, key: str, value: Any) -> Any: + if key in _BOOL_PARAMS: + # accept a real bool, numpy bool, or an integer 0/1; store a Python bool + if isinstance(value, (bool, np.bool_)): + return bool(value) + if isinstance(value, (int, np.integer)): + if int(value) not in (0, 1): + raise ValueError(f"Parameter '{key}' must be 0 or 1 when given as an int.") + return bool(value) + raise TypeError(f"Parameter '{key}' must be a bool (or 0/1).") + + if key in _INT_PARAMS: + # accept any Python/numpy integer or an integer-valued float; store a Python int + if isinstance(value, (bool, np.bool_)): + raise TypeError(f"Parameter '{key}' must be an int.") + if isinstance(value, (int, np.integer)): + value = int(value) + elif isinstance(value, (float, np.floating)) and float(value).is_integer(): + value = int(value) + else: + raise TypeError(f"Parameter '{key}' must be an int.") + if key in _POSITIVE_INT_PARAMS and value <= 0: + raise ValueError(f"Parameter '{key}' must be positive.") + if key not in _POSITIVE_INT_PARAMS and value < 0: + raise ValueError(f"Parameter '{key}' must be nonnegative.") + if value > _INT32_MAX: + raise ValueError(f"Parameter '{key}' must not exceed {_INT32_MAX} (int32 range).") + return value + + if key in _FLOAT_PARAMS: + # accept any real number (Python/numpy int or float); store a Python float + if isinstance(value, (bool, np.bool_)): + raise TypeError(f"Parameter '{key}' must be a number.") + if not isinstance(value, (int, float, np.integer, np.floating)): + raise TypeError(f"Parameter '{key}' must be a number.") + value = float(value) + if not np.isfinite(value): + raise ValueError(f"Parameter '{key}' must be finite.") + if key in _POSITIVE_FLOAT_PARAMS and value <= 0.0: + raise ValueError(f"Parameter '{key}' must be positive.") + if key in _NONNEGATIVE_FLOAT_PARAMS and value < 0.0: + raise ValueError(f"Parameter '{key}' must be nonnegative.") + return value + + if key in _STRING_PARAMS: + if not isinstance(value, str): + raise TypeError(f"Parameter '{key}' must be str.") + value = value.lower() + if value not in ("l2", "linf"): + raise ValueError("Parameter 'optimality_norm' must be 'l2' or 'linf'.") + return value + + return value + def setParam(self, name: str, value: Any) -> None: """ Set the value of a solver parameter by name. """ - key = PDLP._PARAM_ALIAS.get(name, name) - self._params[key] = value + # resolve name and store + key = self._resolve_param_key(name) + self._params[key] = self._validate_param_value(key, value) def getParam(self, name: str) -> Any: """ Get the value of a solver parameter by name. """ - key = PDLP._PARAM_ALIAS.get(name, name) - return self._params.get(key) + # resolve name and return + key = self._resolve_param_key(name) + return self._params[key] def setParams(self, /, **kwargs) -> None: """ Set multiple solver parameters by name. """ + updates = {} for k, v in kwargs.items(): - self.setParam(k, v) + key = self._resolve_param_key(k) + updates[key] = self._validate_param_value(key, v) + self._params.update(updates) + + def resetParams(self) -> None: + """ + Reset all solver parameters to their backend default values. + """ + self._params = dict(self._default_params) def optimize(self): """ @@ -386,6 +612,7 @@ def optimize(self): # check model sense if self.ModelSense not in (PDLP.MINIMIZE, PDLP.MAXIMIZE): raise ValueError("model_sense must be PDLP.MINIMIZE or PDLP.MAXIMIZE") + self._validate_bounds() minimize = self.ModelSense == PDLP.MINIMIZE # call the core solver info = solve_once( @@ -402,20 +629,24 @@ def optimize(self): minimize=minimize, ) # solutions - self._x = np.asarray(info.get("X")) if info.get("X") is not None else None - self._y = np.asarray(info.get("Pi")) if info.get("Pi") is not None else None - self._rc = np.asarray(info.get("RC")) if info.get("RC") is not None else None + x = info.get("X") + y = info.get("Pi") + rc = info.get("RC") + self._x = _readonly_array(np.asarray(x)) if x is not None else None + self._y = _readonly_array(np.asarray(y)) if y is not None else None + self._rc = _readonly_array(np.asarray(rc)) if rc is not None else None # objectives & gaps - primal_obj_eff = info.get("PrimalObj") - dual_obj_eff = info.get("DualObj") - self._objval = primal_obj_eff if primal_obj_eff is not None else None - self._dualobj = dual_obj_eff if dual_obj_eff is not None else None + self._objval = info.get("PrimalObj") + self._dualobj = info.get("DualObj") self._gap = info.get("ObjectiveGap") self._rel_gap = info.get("RelativeObjectiveGap") # status & counters - self._status = str(info.get("Status")) if info.get("Status") is not None else None - self._status_code = int(info.get("StatusCode")) if info.get("StatusCode") is not None else None - self._iter = int(info.get("Iterations")) if info.get("Iterations") is not None else None + status = info.get("Status") + status_code = info.get("StatusCode") + iters = info.get("Iterations") + self._status_name = str(status) if status is not None else None + self._status_code = int(status_code) if status_code is not None else None + self._iter = int(iters) if iters is not None else None self._runtime = info.get("RuntimeSec") self._rescale_time = info.get("RescalingTimeSec") # residuals @@ -424,10 +655,9 @@ def optimize(self): # rays self._max_p_ray = info.get("MaxPrimalRayInfeas") self._max_d_ray = info.get("MaxDualRayInfeas") - p_ray_lin_eff = info.get("PrimalRayLinObj") - d_ray_obj_eff = info.get("DualRayObj") - self._p_ray_lin_obj = p_ray_lin_eff if p_ray_lin_eff is not None else None - self._d_ray_obj = d_ray_obj_eff if d_ray_obj_eff is not None else None + self._p_ray_lin_obj = info.get("PrimalRayLinObj") + self._d_ray_obj = info.get("DualRayObj") + return self def _clear_solution_cache(self) -> None: """ @@ -436,7 +666,7 @@ def _clear_solution_cache(self) -> None: self._x = self._y = self._rc = None self._objval = self._dualobj = None self._gap = self._rel_gap = None - self._status = None + self._status_name = None self._status_code = None self._iter = None self._runtime = self._rescale_time = None @@ -445,6 +675,88 @@ def _clear_solution_cache(self) -> None: self._max_p_ray = self._max_d_ray = None self._p_ray_lin_obj = self._d_ray_obj = None + def _validate_bounds(self) -> None: + _check_bounds(self._lb, self._ub, "variable bounds") + _check_bounds(self._constr_lb, self._constr_ub, "constraint bounds") + + # model data (read/write; assignment reroutes through the validating setters) + @property + def c(self) -> Optional[np.ndarray]: + """Objective coefficient vector.""" + return self._c + + @c.setter + def c(self, value: ArrayLike) -> None: + self.setObjectiveVector(value) + + @property + def c0(self) -> float: + """Objective constant term.""" + return self._c0 + + @c0.setter + def c0(self, value: float) -> None: + self.setObjectiveConstant(value) + + @property + def A(self): + """Constraint matrix (CSR for sparse input, dense ndarray otherwise).""" + return self._A + + @A.setter + def A(self, value) -> None: + self.setConstraintMatrix(value) + + @property + def lb(self) -> Optional[np.ndarray]: + """Variable lower bounds (None means -inf).""" + return self._lb + + @lb.setter + def lb(self, value: Optional[ArrayLike]) -> None: + self.setVariableLowerBound(value) + + @property + def ub(self) -> Optional[np.ndarray]: + """Variable upper bounds (None means +inf).""" + return self._ub + + @ub.setter + def ub(self, value: Optional[ArrayLike]) -> None: + self.setVariableUpperBound(value) + + @property + def constr_lb(self) -> Optional[np.ndarray]: + """Constraint lower bounds (None means -inf).""" + return self._constr_lb + + @constr_lb.setter + def constr_lb(self, value: Optional[ArrayLike]) -> None: + self.setConstraintLowerBound(value) + + @property + def constr_ub(self) -> Optional[np.ndarray]: + """Constraint upper bounds (None means +inf).""" + return self._constr_ub + + @constr_ub.setter + def constr_ub(self, value: Optional[ArrayLike]) -> None: + self.setConstraintUpperBound(value) + + @property + def ModelSense(self) -> int: + """Objective sense: PDLP.MINIMIZE or PDLP.MAXIMIZE.""" + return self._model_sense + + @ModelSense.setter + def ModelSense(self, value: int) -> None: + # validate sense + if value not in (PDLP.MINIMIZE, PDLP.MAXIMIZE): + raise ValueError("ModelSense must be PDLP.MINIMIZE or PDLP.MAXIMIZE") + self._model_sense = value + # clear cached solution + self._clear_solution_cache() + @property def X(self) -> Optional[np.ndarray]: return self._x @@ -474,12 +786,17 @@ def RelGap(self) -> Optional[float]: return self._rel_gap @property - def Status(self) -> Optional[str]: - return self._status + def Status(self) -> Optional[int]: + """ + Integer termination status code. Compare against the constants in + cupdlpx.PDLP, e.g. ``model.Status == PDLP.OPTIMAL``. + """ + return self._status_code @property - def StatusCode(self) -> Optional[int]: - return self._status_code + def StatusName(self) -> Optional[str]: + """Human-readable termination status name, e.g. ``'OPTIMAL'``.""" + return self._status_name @property def IterCount(self) -> Optional[int]: @@ -519,8 +836,10 @@ def DualRayObj(self) -> Optional[float]: @property def PrimalInfeas(self) -> Optional[float]: + """Alias of RelPrimalResidual (relative, not absolute infeasibility).""" return self._rel_p_res @property def DualInfeas(self) -> Optional[float]: - return self._rel_d_res \ No newline at end of file + """Alias of RelDualResidual (relative, not absolute infeasibility).""" + return self._rel_d_res diff --git a/python_bindings/_core_bindings.cpp b/python_bindings/_core_bindings.cpp index f2919de..40971e7 100644 --- a/python_bindings/_core_bindings.cpp +++ b/python_bindings/_core_bindings.cpp @@ -15,6 +15,7 @@ limitations under the License. */ #include +#include #include #include #include @@ -226,12 +227,38 @@ static int status_to_code(termination_reason_t r) return 4; case TERMINATION_REASON_INFEASIBLE_OR_UNBOUNDED: return 5; + case TERMINATION_REASON_FEAS_POLISH_SUCCESS: + return 6; case TERMINATION_REASON_UNSPECIFIED: default: return -1; } } +static void validate_result_dimensions(const cupdlpx_result_t *res, int expected_n, int expected_m) +{ + if (res->num_variables != expected_n || res->num_constraints != expected_m) + { + throw std::runtime_error("solve_lp_problem returned result dimensions " + + std::to_string(res->num_variables) + "x" + + std::to_string(res->num_constraints) + + ", expected " + std::to_string(expected_n) + + "x" + std::to_string(expected_m) + "."); + } + if (expected_n > 0 && !res->primal_solution) + { + throw std::runtime_error("solve_lp_problem returned NULL primal_solution."); + } + if (expected_m > 0 && !res->dual_solution) + { + throw std::runtime_error("solve_lp_problem returned NULL dual_solution."); + } + if (expected_n > 0 && !res->reduced_cost) + { + throw std::runtime_error("solve_lp_problem returned NULL reduced_cost."); + } +} + // get default parameters as Python dict static py::dict get_default_params_py() { @@ -294,17 +321,42 @@ static void parse_params_from_python(py::object params_obj, pdhg_parameters_t *p auto getf = [&](const char *k, double &tgt) { if (d.contains(k)) - tgt = py::cast(d[k]); + { + py::object val = d[k]; + if (py::isinstance(val)) + { + throw std::invalid_argument(std::string(k) + " must be a number."); + } + tgt = py::cast(val); + if (!std::isfinite(tgt)) + { + throw std::invalid_argument(std::string(k) + " must be finite."); + } + } }; auto geti = [&](const char *k, int &tgt) { if (d.contains(k)) - tgt = py::cast(d[k]); + { + py::object val = d[k]; + if (py::isinstance(val) || !py::isinstance(val)) + { + throw std::invalid_argument(std::string(k) + " must be an int."); + } + tgt = py::cast(val); + } }; auto getb = [&](const char *k, bool &tgt) { if (d.contains(k)) - tgt = py::cast(d[k]); + { + py::object val = d[k]; + if (!py::isinstance(val)) + { + throw std::invalid_argument(std::string(k) + " must be a bool."); + } + tgt = py::cast(val); + } }; auto get_norm = [&](const char *k, norm_type_t &tgt) { @@ -364,6 +416,131 @@ static void parse_params_from_python(py::object params_obj, pdhg_parameters_t *p getb("presolve", p->presolve); getf("matrix_zero_tol", p->matrix_zero_tol); + + if (p->termination_evaluation_frequency <= 0) + throw std::invalid_argument("termination_evaluation_frequency must be positive."); + if (p->termination_criteria.iteration_limit < 0) + throw std::invalid_argument("iteration_limit must be nonnegative."); + if (p->l_inf_ruiz_iterations < 0) + throw std::invalid_argument("l_inf_ruiz_iterations must be nonnegative."); + if (p->sv_max_iter <= 0) + throw std::invalid_argument("sv_max_iter must be positive."); + if (p->termination_criteria.eps_optimal_relative <= 0.0) + throw std::invalid_argument("eps_optimal_relative must be positive."); + if (p->termination_criteria.eps_feasible_relative <= 0.0) + throw std::invalid_argument("eps_feasible_relative must be positive."); + if (p->termination_criteria.eps_feas_polish_relative <= 0.0) + throw std::invalid_argument("eps_feas_polish_relative must be positive."); + if (p->sv_tol <= 0.0) + throw std::invalid_argument("sv_tol must be positive."); + if (p->termination_criteria.time_sec_limit < 0.0) + throw std::invalid_argument("time_sec_limit must be nonnegative."); + if (p->matrix_zero_tol < 0.0) + throw std::invalid_argument("matrix_zero_tol must be nonnegative."); +} + +// throw if a 1D array's length differs from the expected value +static void expect_len(py::object obj, py::ssize_t expected, const char *name) +{ + py::array arr = py::cast(obj); + if (arr.ndim() != 1) + { + throw std::invalid_argument(std::string(name) + " must be 1D."); + } + if (arr.size() != expected) + { + throw std::invalid_argument(std::string(name) + " has wrong length: expected " + + std::to_string(expected) + ", got " + std::to_string((long long)arr.size())); + } +} + +static void validate_finite_array(const double *data, py::ssize_t size, const char *name) +{ + if (!data) + { + return; + } + for (py::ssize_t i = 0; i < size; ++i) + { + if (!std::isfinite(data[i])) + { + throw std::invalid_argument(std::string(name) + " must contain only finite values."); + } + } +} + +static void validate_no_nan_array(const double *data, py::ssize_t size, const char *name) +{ + if (!data) + { + return; + } + for (py::ssize_t i = 0; i < size; ++i) + { + if (std::isnan(data[i])) + { + throw std::invalid_argument(std::string(name) + " must not contain NaN."); + } + } +} + +static void validate_bounds(const double *lower, const double *upper, int size, const char *name) +{ + if (!lower || !upper) + { + return; + } + for (int i = 0; i < size; ++i) + { + if (lower[i] > upper[i]) + { + throw std::invalid_argument(std::string(name) + ": lower bounds must be <= upper bounds."); + } + } +} + +// validate a compressed (CSR/CSC) index structure +static void validate_compressed(const int32_t *indptr, const int32_t *indices, int major, int minor, int nnz, + const char *fmt) +{ + if (indptr[0] != 0) + { + throw std::invalid_argument(std::string(fmt) + ".indptr[0] must be 0."); + } + for (int i = 0; i < major; ++i) + { + if (indptr[i] > indptr[i + 1]) + { + throw std::invalid_argument(std::string(fmt) + ".indptr must be non-decreasing."); + } + } + if (indptr[major] != nnz) + { + throw std::invalid_argument(std::string(fmt) + ".indptr[-1] must equal nnz."); + } + for (int k = 0; k < nnz; ++k) + { + if (indices[k] < 0 || indices[k] >= minor) + { + throw std::invalid_argument(std::string(fmt) + " has an index out of range [0, dim)."); + } + } +} + +// validate COO row/column indices are within [0, m) and [0, n) respectively +static void validate_coo(const int32_t *row, const int32_t *col, int m, int n, int nnz) +{ + for (int k = 0; k < nnz; ++k) + { + if (row[k] < 0 || row[k] >= m) + { + throw std::invalid_argument("coo.row has an index out of range [0, m)."); + } + if (col[k] < 0 || col[k] >= n) + { + throw std::invalid_argument("coo.col has an index out of range [0, n)."); + } + } } // view of matrix from Python @@ -394,6 +571,7 @@ static PyMatrixView get_matrix_from_python(py::object A) { throw std::invalid_argument("dense matrix must be 2D"); } + validate_finite_array(static_cast(req.ptr), d.size(), "dense matrix"); desc.m = static_cast(req.shape[0]); desc.n = static_cast(req.shape[1]); desc.fmt = matrix_dense; @@ -413,11 +591,17 @@ static PyMatrixView get_matrix_from_python(py::object A) py::object ci = A.attr("indices"); py::object vv = A.attr("data"); py::array v64 = get_array_f64_c_contig(vv, "csr.data(float64)"); // get contiguous data array + validate_finite_array(static_cast(v64.request().ptr), v64.size(), "csr.data"); desc.fmt = matrix_csr; desc.data.csr.nnz = static_cast(v64.size()); + // check index array lengths before dereferencing + expect_len(rp, static_cast(desc.m) + 1, "csr.indptr"); + expect_len(ci, static_cast(desc.data.csr.nnz), "csr.indices"); desc.data.csr.row_ptr = get_index_ptr_i32(rp, "csr.indptr", out.keep, out.keep.tmp_rowptr); desc.data.csr.col_ind = get_index_ptr_i32(ci, "csr.indices", out.keep, out.keep.tmp_colind); desc.data.csr.vals = static_cast(v64.request().ptr); + // validate structure (protects the public solve_once entry) + validate_compressed(desc.data.csr.row_ptr, desc.data.csr.col_ind, desc.m, desc.n, desc.data.csr.nnz, "csr"); out.keep.owners.push_back(v64); // keep alive return out; } @@ -428,11 +612,17 @@ static PyMatrixView get_matrix_from_python(py::object A) py::object ri = A.attr("indices"); py::object vv = A.attr("data"); py::array v64 = get_array_f64_c_contig(vv, "csc.data(float64)"); // get contiguous data array + validate_finite_array(static_cast(v64.request().ptr), v64.size(), "csc.data"); desc.fmt = matrix_csc; desc.data.csc.nnz = static_cast(v64.size()); + // check index array lengths before dereferencing + expect_len(cp, static_cast(desc.n) + 1, "csc.indptr"); + expect_len(ri, static_cast(desc.data.csc.nnz), "csc.indices"); desc.data.csc.col_ptr = get_index_ptr_i32(cp, "csc.indptr", out.keep, out.keep.tmp_rowptr); desc.data.csc.row_ind = get_index_ptr_i32(ri, "csc.indices", out.keep, out.keep.tmp_colind); desc.data.csc.vals = static_cast(v64.request().ptr); + // validate structure (major=n, minor=m for CSC) + validate_compressed(desc.data.csc.col_ptr, desc.data.csc.row_ind, desc.n, desc.m, desc.data.csc.nnz, "csc"); out.keep.owners.push_back(v64); // keep alive return out; } @@ -443,11 +633,17 @@ static PyMatrixView get_matrix_from_python(py::object A) py::object cc = A.attr("col"); py::object vv = A.attr("data"); py::array v64 = get_array_f64_c_contig(vv, "coo.data(float64)"); // get contiguous data array + validate_finite_array(static_cast(v64.request().ptr), v64.size(), "coo.data"); desc.fmt = matrix_coo; desc.data.coo.nnz = static_cast(v64.size()); + // check index array lengths before dereferencing + expect_len(rr, static_cast(desc.data.coo.nnz), "coo.row"); + expect_len(cc, static_cast(desc.data.coo.nnz), "coo.col"); desc.data.coo.row_ind = get_index_ptr_i32(rr, "coo.row", out.keep, out.keep.tmp_row); desc.data.coo.col_ind = get_index_ptr_i32(cc, "coo.col", out.keep, out.keep.tmp_col); desc.data.coo.vals = static_cast(v64.request().ptr); + // validate indices are within [0, m) x [0, n) + validate_coo(desc.data.coo.row_ind, desc.data.coo.col_ind, desc.m, desc.n, desc.data.coo.nnz); out.keep.owners.push_back(v64); // keep alive return out; } @@ -460,7 +656,7 @@ static PyMatrixView get_matrix_from_python(py::object A) static py::dict solve_once(py::object A, py::object objective_vector, // c py::object objective_constant, // c0 (optional → 0) - py::object variable_lower_bound, // lb (optional → 0) + py::object variable_lower_bound, // lb (optional → -inf) py::object variable_upper_bound, // ub (optional → inf) py::object constraint_lower_bound, // l (optional → -inf) py::object constraint_upper_bound, // u (optional → inf) @@ -485,12 +681,23 @@ static py::dict solve_once(py::object A, const double *ub_ptr = get_arr_ptr_f64_or_null(variable_upper_bound, "variable_upper_bound", view.keep); const double *l_ptr = get_arr_ptr_f64_or_null(constraint_lower_bound, "constraint_lower_bound", view.keep); const double *u_ptr = get_arr_ptr_f64_or_null(constraint_upper_bound, "constraint_upper_bound", view.keep); + validate_finite_array(c_ptr, n, "objective_vector"); + validate_no_nan_array(lb_ptr, n, "variable_lower_bound"); + validate_no_nan_array(ub_ptr, n, "variable_upper_bound"); + validate_no_nan_array(l_ptr, m, "constraint_lower_bound"); + validate_no_nan_array(u_ptr, m, "constraint_upper_bound"); + validate_bounds(lb_ptr, ub_ptr, n, "variable bounds"); + validate_bounds(l_ptr, u_ptr, m, "constraint bounds"); // get objective constant double c0_local = 0.0; double *c0_ptr = nullptr; if (objective_constant && !objective_constant.is_none()) { c0_local = py::cast(objective_constant); + if (!std::isfinite(c0_local)) + { + throw std::invalid_argument("objective_constant must be finite."); + } c0_ptr = &c0_local; } @@ -509,6 +716,8 @@ static py::dict solve_once(py::object A, { throw std::runtime_error("create_lp_problem failed."); } + // free the problem even if a conversion/validation below throws + std::unique_ptr prob_guard(prob, &lp_problem_free); // set warm start values if provided if ((primal_start && !primal_start.is_none()) || (dual_start && !dual_start.is_none())) @@ -518,6 +727,8 @@ static py::dict solve_once(py::object A, ensure_len_or_null(dual_start, "dual_start", m); const double *primal_ptr = get_arr_ptr_f64_or_null(primal_start, "primal_start", view.keep); const double *dual_ptr = get_arr_ptr_f64_or_null(dual_start, "dual_start", view.keep); + validate_finite_array(primal_ptr, n, "primal_start"); + validate_finite_array(dual_ptr, m, "dual_start"); set_start_values(prob, primal_ptr, dual_ptr); } @@ -532,13 +743,17 @@ static py::dict solve_once(py::object A, py::gil_scoped_release release; res = solve_lp_problem(prob, &local_params); } - lp_problem_free(prob); + // problem is no longer needed once the solve returns + prob_guard.reset(); if (!res) { throw std::runtime_error("solve_lp_problem returned NULL."); } + // free the result even if a conversion below throws + std::unique_ptr res_guard(res, &cupdlpx_result_free); // parse result + validate_result_dimensions(res, n, m); const int n_out = res->num_variables; const int m_out = res->num_constraints; py::array_t x({n_out}); @@ -546,9 +761,15 @@ static py::dict solve_once(py::object A, py::array_t rc({n_out}); { auto xb = x.request(), yb = y.request(), rcb = rc.request(); - std::memcpy(xb.ptr, res->primal_solution, sizeof(double) * n_out); - std::memcpy(yb.ptr, res->dual_solution, sizeof(double) * m_out); - std::memcpy(rcb.ptr, res->reduced_cost, sizeof(double) * n_out); + if (n_out > 0) + { + std::memcpy(xb.ptr, res->primal_solution, sizeof(double) * n_out); + std::memcpy(rcb.ptr, res->reduced_cost, sizeof(double) * n_out); + } + if (m_out > 0) + { + std::memcpy(yb.ptr, res->dual_solution, sizeof(double) * m_out); + } } // build info dict py::dict info; @@ -576,9 +797,7 @@ static py::dict solve_once(py::object A, info["PrimalRayLinObj"] = res->primal_ray_linear_objective; info["DualRayObj"] = res->dual_ray_objective; - // free result - cupdlpx_result_free(res); - + // res freed by res_guard on return return info; } diff --git a/test/test_api_surface.py b/test/test_api_surface.py new file mode 100644 index 0000000..b0e19e3 --- /dev/null +++ b/test/test_api_surface.py @@ -0,0 +1,429 @@ +# Copyright 2025 Haihao Lu +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +import os + +import numpy as np +import scipy.sparse as sp +import pytest + +import cupdlpx +from cupdlpx import Model, PDLP, read +from cupdlpx._core import solve_once + + +def test_public_exports_include_documented_api(): + assert cupdlpx.PDLP is PDLP + assert cupdlpx.Model is Model + assert cupdlpx.read is read + assert isinstance(cupdlpx.__version__, str) + assert set(cupdlpx.__all__) == {"Model", "PDLP", "read", "__version__"} + + +def _model(base_lp_data): + c, A, l, u, lb, ub = base_lp_data + model = Model(c, A, l, u, lb, ub) + model.setParams(OutputFlag=False, Presolve=False) + return model + + +def test_data_property_reads(base_lp_data): + """Read-side of the model-data properties.""" + model = _model(base_lp_data) + assert model.c.shape == (2,) + assert model.c0 == 0.0 + assert model.A.shape == (3, 2) + assert model.constr_lb.shape == (3,) + assert model.constr_ub.shape == (3,) + assert model.lb is None + assert model.ub is None + assert model.ModelSense == PDLP.MINIMIZE + + +def test_data_property_writes(base_lp_data): + """Direct assignment routes through the validating setters.""" + model = _model(base_lp_data) + model.c = [2.0, 3.0] + assert np.allclose(model.c, [2.0, 3.0]) + model.c0 = 4.0 + assert model.c0 == 4.0 + model.lb = [0.0, 0.0] + assert np.allclose(model.lb, [0.0, 0.0]) + model.ub = [10.0, 10.0] + assert np.allclose(model.ub, [10.0, 10.0]) + model.constr_lb = [1.0, 2.0, 3.0] + assert np.allclose(model.constr_lb, [1.0, 2.0, 3.0]) + model.constr_ub = [4.0, 5.0, 6.0] + assert np.allclose(model.constr_ub, [4.0, 5.0, 6.0]) + model.A = sp.csr_matrix(np.array([[1.0, 0.0], [0.0, 1.0], [1.0, 1.0]])) + assert model.A.shape == (3, 2) + model.ModelSense = PDLP.MAXIMIZE + assert model.ModelSense == PDLP.MAXIMIZE + + +def test_bounds_none_clears(base_lp_data): + """Passing None clears optional bounds.""" + model = _model(base_lp_data) + model.lb = [0.0, 0.0] + model.lb = None + assert model.lb is None + model.ub = None + assert model.ub is None + model.constr_lb = None + assert model.constr_lb is None + model.constr_ub = None + assert model.constr_ub is None + + +def test_property_validation_errors(base_lp_data): + """Invalid assignments raise immediately instead of failing in the solver.""" + model = _model(base_lp_data) + with pytest.raises(ValueError): + model.c = [1.0, 2.0, 3.0] # wrong length + with pytest.raises(ValueError): + model.c = [[1.0, 2.0]] # not 1D + with pytest.raises(TypeError): + model.A = "not a matrix" + with pytest.raises(ValueError): + model.A = np.ones(3) # not 2D + with pytest.raises(ValueError): + model.A = np.ones((3, 5)) # wrong number of variables + with pytest.raises(ValueError): + model.lb = [0.0] # wrong length + with pytest.raises(ValueError): + model.ub = [0.0] + with pytest.raises(ValueError): + model.constr_lb = [0.0, 0.0] # need 3 + with pytest.raises(ValueError): + model.constr_ub = [0.0, 0.0] + with pytest.raises(ValueError): + model.ModelSense = 99 + + +def test_bound_nan_validation(base_lp_data): + model = _model(base_lp_data) + with pytest.raises(ValueError): + model.lb = [0.0, np.nan] + with pytest.raises(ValueError): + model.ub = [np.nan, 1.0] + with pytest.raises(ValueError): + model.constr_lb = [0.0, np.nan, -np.inf] + with pytest.raises(ValueError): + model.constr_ub = [0.0, 1.0, np.nan] + + model.lb = [0.0, -np.inf] + model.ub = [np.inf, 1.0] + model.constr_lb = [5.0, -np.inf, -np.inf] + model.constr_ub = [5.0, np.inf, 8.0] + + +def test_invalid_objective_assignment_preserves_previous_value(base_lp_data): + model = _model(base_lp_data) + original = model.c.copy() + with pytest.raises(ValueError): + model.c = [1.0, 2.0, 3.0] + assert np.allclose(model.c, original) + + +def test_dense_inputs_are_copied(base_lp_data): + c, A, l, u, lb, ub = base_lp_data + c = c.copy() + A = A.copy() + l = l.copy() + u = u.copy() + model = Model(c, A, l, u, lb, ub) + c[0] = 99.0 + A[0, 0] = 99.0 + l[0] = 99.0 + u[0] = 99.0 + assert np.allclose(model.c, [1.0, 1.0]) + assert model.A[0, 0] == 1.0 + assert model.constr_lb[0] == 5.0 + assert model.constr_ub[0] == 5.0 + + +def test_model_arrays_are_read_only(base_lp_data): + c, A, l, u, lb, ub = base_lp_data + model = Model(c, A, l, u, lb, ub) + with pytest.raises(ValueError): + model.c[0] = 99.0 + with pytest.raises(ValueError): + model.A[0, 0] = 99.0 + with pytest.raises(ValueError): + model.constr_lb[0] = 99.0 + with pytest.raises(ValueError): + model.constr_ub[0] = 99.0 + + model.c = [2.0, 3.0] + assert np.allclose(model.c, [2.0, 3.0]) + + +def test_sparse_model_arrays_are_read_only(base_lp_data): + c, A, l, u, lb, ub = base_lp_data + model = Model(c, sp.csr_matrix(A), l, u, lb, ub) + with pytest.raises(ValueError): + model.A.data[0] = 99.0 + with pytest.raises(ValueError): + model.A.indices[0] = 99 + with pytest.raises(ValueError): + model.A.indptr[0] = 99 + + +def test_bound_order_validation(base_lp_data): + model = _model(base_lp_data) + model.lb = [2.0, 0.0] + model.ub = [1.0, 1.0] + with pytest.raises(ValueError): + model.optimize() + + with pytest.raises(ValueError): + Model( + [1.0], + np.array([[1.0]]), + constraint_lower_bound=[2.0], + constraint_upper_bound=[1.0], + ) + + +def test_set_params_is_transactional(base_lp_data): + model = _model(base_lp_data) + old_time_limit = model.getParam("TimeLimit") + with pytest.raises(KeyError): + model.setParams(TimeLimit=123.0, DefinitelyNotAParam=1) + assert model.getParam("TimeLimit") == old_time_limit + + +def test_param_value_validation(base_lp_data): + model = _model(base_lp_data) + # float params accept ints and numpy numbers + model.setParam("TimeLimit", 12) + assert model.getParam("TimeLimit") == 12.0 + model.setParam("TimeLimit", np.float64(30.0)) + assert model.getParam("TimeLimit") == 30.0 + model.setParam("OptimalityNorm", "LINF") + assert model.getParam("OptimalityNorm") == "linf" + + # bool params accept real/numpy bools and 0/1 (coerced to a Python bool) + model.setParam("OutputFlag", 1) + assert model.getParam("OutputFlag") is True + model.setParam("OutputFlag", 0) + assert model.getParam("OutputFlag") is False + model.setParam("OutputFlag", np.bool_(True)) + assert model.getParam("OutputFlag") is True + + # int params accept numpy ints and integer-valued floats (coerced to a Python int) + model.setParam("IterationLimit", np.int64(1000)) + assert model.getParam("IterationLimit") == 1000 and isinstance(model.getParam("IterationLimit"), int) + model.setParam("IterationLimit", 2000.0) + assert model.getParam("IterationLimit") == 2000 + + # still-invalid inputs + with pytest.raises(TypeError): + model.setParam("IterationLimit", False) # bool is not an int here + with pytest.raises(TypeError): + model.setParam("IterationLimit", 1.5) # non-integer float + with pytest.raises(TypeError): + model.setParam("OutputFlag", "yes") # string is not a bool + with pytest.raises(ValueError): + model.setParam("OutputFlag", 2) # only 0/1 allowed as int + with pytest.raises(ValueError): + model.setParam("IterationLimit", -1) + with pytest.raises(ValueError): + model.setParam("TermCheckFreq", 0) + with pytest.raises(ValueError): + model.setParam("OptimalityTol", 0.0) + with pytest.raises(ValueError): + model.setParam("TimeLimit", -1.0) + with pytest.raises(ValueError): + model.setParam("OptimalityNorm", "l1") + + +def test_direct_core_param_value_validation(base_lp_data): + c, A, l, u, lb, ub = base_lp_data + with pytest.raises(ValueError): + solve_once(A, c, None, lb, ub, l, u, params={"verbose": 1}) + with pytest.raises(ValueError): + solve_once(A, c, None, lb, ub, l, u, params={"iteration_limit": False}) + with pytest.raises(ValueError): + solve_once(A, c, None, lb, ub, l, u, params={"eps_optimal_relative": 0.0}) + with pytest.raises(ValueError): + solve_once(A, c, None, lb, ub, l, u, params={"time_sec_limit": float("nan")}) + with pytest.raises(ValueError): + solve_once(A, c, None, lb, ub, l, u, params={"optimality_norm": "l1"}) + + +def test_direct_core_model_data_validation(base_lp_data): + c, A, l, u, lb, ub = base_lp_data + with pytest.raises(ValueError): + solve_once(A, [np.nan, 1.0], None, lb, ub, l, u) + with pytest.raises(ValueError): + solve_once([[np.nan, 1.0], [0.0, 1.0], [3.0, 2.0]], c, None, lb, ub, l, u) + with pytest.raises(ValueError): + solve_once(A, c, np.nan, lb, ub, l, u) + with pytest.raises(ValueError): + solve_once(A, c, None, [0.0, np.nan], ub, l, u) + with pytest.raises(ValueError): + solve_once(A, c, None, [2.0, 0.0], [1.0, 1.0], l, u) + with pytest.raises(ValueError): + solve_once(A, c, None, lb, ub, l, u, primal_start=[np.nan, 0.0]) + + +def test_set_params_value_validation_is_transactional(base_lp_data): + model = _model(base_lp_data) + old_time_limit = model.getParam("TimeLimit") + with pytest.raises(TypeError): + model.setParams(TimeLimit=123.0, OutputFlag="yes") # OutputFlag invalid + assert model.getParam("TimeLimit") == old_time_limit + + +def test_reset_params(base_lp_data): + c, A, l, u, lb, ub = base_lp_data + model = Model(c, A, l, u, lb, ub) + default_time_limit = model.getParam("TimeLimit") + default_output_flag = model.getParam("OutputFlag") + + model.setParams(TimeLimit=12.0, OutputFlag=not default_output_flag) + model.resetParams() + assert model.getParam("TimeLimit") == default_time_limit + assert model.getParam("OutputFlag") is default_output_flag + + +def test_matrix_row_change_conflicts_with_bounds(base_lp_data): + """Reassigning A with a different row count than existing bounds raises.""" + model = _model(base_lp_data) + # existing constraint bounds have 3 rows; a 2-row matrix must be rejected + # (constraint_lower_bound check) + two_row = sp.csr_matrix(np.array([[1.0, 0.0], [0.0, 1.0]])) + with pytest.raises(ValueError): + model.A = two_row + # now clear the lower bound so the upper-bound check is the one that fires + model.constr_lb = None + with pytest.raises(ValueError): + model.A = two_row + + +def test_params_view(base_lp_data): + """_ParamsView get/set via attribute and mapping-style access.""" + model = _model(base_lp_data) + model.Params.TimeLimit = 123.0 + assert model.Params.TimeLimit == 123.0 + model.Params["OptimalityTol"] = 1e-5 + assert model.Params["OptimalityTol"] == 1e-5 + assert model.getParam("TimeLimit") == 123.0 + assert "TimeLimit" in model.Params + assert "time_sec_limit" in model.Params + assert "DefinitelyNotAParam" not in model.Params + assert len(list(model.Params.keys())) > 0 + assert len(list(model.Params.values())) == len(list(model.Params.keys())) + assert dict(model.Params.items())["time_sec_limit"] == 123.0 + with pytest.raises(AttributeError): + _ = model.Params.NoSuchParameter + + +def test_param_aliases_resolve_consistently(base_lp_data): + """Every public parameter alias resolves to the same backend parameter.""" + model = _model(base_lp_data) + for alias, key in PDLP._PARAM_ALIAS.items(): + assert alias in model.Params + assert key in model.Params + value = model.getParam(key) + assert model.getParam(alias) == value + model.setParam(alias, value) + assert getattr(model.Params, alias) == value + assert model.Params[alias] == value + + +def test_param_name_validation(base_lp_data): + """Unknown parameter names are rejected in both set and get.""" + model = _model(base_lp_data) + with pytest.raises(KeyError): + model.setParam("DefinitelyNotAParam", 1) + with pytest.raises(KeyError): + model.getParam("DefinitelyNotAParam") + + +def test_non_contiguous_and_typed_inputs(base_lp_data): + """Helper conversions handle non-contiguous dense and non-float64 sparse.""" + _, _, l, u, _, _ = base_lp_data + # non-C-contiguous objective vector (a strided column view) + c_nc = np.ones((2, 4))[:, 1] + assert not c_nc.flags["C_CONTIGUOUS"] + # integer-typed sparse matrix (exercises the float64 upcast path) + A = sp.csr_matrix(np.array([[1, 2], [0, 1], [3, 2]], dtype=np.int64)) + model = Model(c_nc, A, l, u) + assert model.c.dtype == np.float64 + assert model.A.dtype == np.float64 + assert model.A.indices.dtype == np.int32 + + +def test_int64_index_downcast(base_lp_data): + """A float64 CSR with int64 indices is downcast to int32 without mutating it.""" + _, _, l, u, _, _ = base_lp_data + c = np.array([1.0, 1.0]) + A = sp.csr_matrix(np.array([[1.0, 2.0], [0.0, 1.0], [3.0, 2.0]])) + A.indices = A.indices.astype(np.int64) + A.indptr = A.indptr.astype(np.int64) + model = Model(c, A, l, u) + assert model.A.indices.dtype == np.int32 + assert model.A.indptr.dtype == np.int32 + # caller's matrix must be untouched + assert A.indices.dtype == np.int64 + assert A.indptr.dtype == np.int64 + + +def test_init_requires_2d_matrix(): + with pytest.raises(ValueError): + Model(np.ones(2), np.ones(3), None, None) # 1D "matrix" + + +def test_solution_attributes_accessible(base_lp_data): + """Every result attribute is readable after optimize().""" + model = _model(base_lp_data) + model.optimize() + # touching each property exercises its getter + attrs = [ + model.X, model.Pi, model.RC, model.ObjVal, model.DualObj, + model.Gap, model.RelGap, model.Status, model.StatusName, + model.IterCount, model.Runtime, model.RescalingTime, + model.RelPrimalResidual, model.RelDualResidual, + model.MaxPrimalRayInfeas, model.MaxDualRayInfeas, + model.PrimalRayLinObj, model.DualRayObj, + model.PrimalInfeas, model.DualInfeas, + ] + assert model.Status == PDLP.OPTIMAL + assert model.StatusName == "OPTIMAL" + assert len(attrs) == 20 + + +def test_optimize_rejects_bad_sense(base_lp_data): + """The defensive sense check inside optimize() rejects a corrupted sense.""" + model = _model(base_lp_data) + model._model_sense = 999 # bypass the property to hit optimize()'s guard + with pytest.raises(ValueError): + model.optimize() + + +def test_read_mps_roundtrip(): + """read() loads an MPS file into a usable Model.""" + path = os.path.join(os.path.dirname(__file__), "cplex2.mps") + model = read(path) + assert model.num_vars > 0 + assert model.num_constrs > 0 + assert model.A.shape == (model.num_constrs, model.num_vars) + assert model.ModelSense in (PDLP.MINIMIZE, PDLP.MAXIMIZE) + + +def test_read_mps_missing_file(): + with pytest.raises(FileNotFoundError): + read("this_file_does_not_exist.mps") diff --git a/test/test_basic.py b/test/test_basic.py index 8ffee5f..a21591e 100644 --- a/test/test_basic.py +++ b/test/test_basic.py @@ -24,8 +24,8 @@ def test_smoke_optimize_runs(base_lp_data): model = Model(c, A, l, u, lb, ub) # turn off output model.setParams(OutputFlag=False, Presolve=False) - # optimize - model.optimize() + # optimize + assert model.optimize() is model def test_minimize_solution_correct(base_lp_data, atol): """ @@ -47,7 +47,7 @@ def test_minimize_solution_correct(base_lp_data, atol): model.optimize() # check status assert hasattr(model, "Status"), "Model.Status not exposed." - assert model.Status == "OPTIMAL", f"Unexpected termination status: {model.Status}" + assert model.Status == PDLP.OPTIMAL, f"Unexpected termination status: {model.Status}" # check primal solution assert hasattr(model, "X"), "Model.X (primal solution) not exposed." assert np.allclose(model.X, [1, 2], atol=atol), f"Unexpected primal solution: {model.X}" @@ -76,17 +76,14 @@ def test_maximize_solution_correct(base_lp_data, atol): c, A, l, u, lb, ub = base_lp_data model = Model(c, A, l, u, lb, ub) # model sense - try: - model.ModelSense = PDLP.MAXIMIZE - except Exception as e: - print(f"cuPDLPx: failed to set model sense to MAXIMIZE.") + model.ModelSense = PDLP.MAXIMIZE # turn off output model.setParams(OutputFlag=False, Presolve=False) # optimize model.optimize() # check status assert hasattr(model, "Status"), "Model.Status not exposed." - assert model.Status == "OPTIMAL", f"Unexpected termination status: {model.Status}" + assert model.Status == PDLP.OPTIMAL, f"Unexpected termination status: {model.Status}" # check primal solution assert hasattr(model, "X"), "Model.X (primal solution) not exposed." assert np.allclose(model.X, [1.5, 1.75], atol=atol), f"Unexpected primal solution: {model.X}" @@ -95,4 +92,4 @@ def test_maximize_solution_correct(base_lp_data, atol): assert np.allclose(model.Pi, [0.25, 0, 0.25], atol=atol), f"Unexpected dual solution: {model.Pi}" # check objective assert hasattr(model, "ObjVal"), "Model.ObjVal (objective value) not exposed." - assert np.isclose(model.ObjVal, 3.25, atol=atol), f"Unexpected objective value: {model.ObjVal}" \ No newline at end of file + assert np.isclose(model.ObjVal, 3.25, atol=atol), f"Unexpected objective value: {model.ObjVal}" diff --git a/test/test_feasibility_polishing.py b/test/test_feasibility_polishing.py index ef63455..6ff7481 100644 --- a/test/test_feasibility_polishing.py +++ b/test/test_feasibility_polishing.py @@ -50,7 +50,7 @@ def test_feasibility_polishing(): m.optimize() # 4. Sanity checks - assert m.Status == "OPTIMAL", f"unexpected status: {m.Status}" + assert m.Status == PDLP.OPTIMAL, f"unexpected status: {m.Status}" assert hasattr(m, "X") and hasattr(m, "ObjVal") # 5. Feasibility-quality check: max violation < 1e-10 diff --git a/test/test_infeasible_unbounded.py b/test/test_infeasible_unbounded.py index 77f01e6..8fffaab 100644 --- a/test/test_infeasible_unbounded.py +++ b/test/test_infeasible_unbounded.py @@ -39,8 +39,7 @@ def test_infeasible_lp(base_lp_data, atol): #model.optimize() # check status #assert hasattr(model, "Status"), "Model.Status not exposed." - #assert model.Status == "PRIMAL_INFEASIBLE", f"Unexpected termination status: {model.Status}" - #assert model.StatusCode == PDLP.PRIMAL_INFEASIBLE, f"Unexpected termination status code: {model.StatusCode}" + #assert model.Status == PDLP.PRIMAL_INFEASIBLE, f"Unexpected termination status: {model.Status}" # check dual ray #assert model.DualRayObj > atol, f"DualRayObj should be positive for dual infeasible, got {model.DualRayObj}" @@ -65,13 +64,10 @@ def test_unbounded_lp(base_lp_data, atol): model.setVariableLowerBound(lb) # turn off output model.setParams(OutputFlag=False, Presolve=False) - # set infeasible tolerance - model.setParams(InfeasibleTol=1e-6) # optimize #model.optimize() # check status #assert hasattr(model, "Status"), "Model.Status not exposed." - #assert model.Status == "DUAL_INFEASIBLE", f"Unexpected termination status: {model.Status}" - #assert model.StatusCode == PDLP.DUAL_INFEASIBLE, f"Unexpected termination status code: {model.StatusCode}" + #assert model.Status == PDLP.DUAL_INFEASIBLE, f"Unexpected termination status: {model.Status}" # check primal ray #assert model.PrimalRayLinObj < -atol, f"PrimalRayLinObj should be negative for dual infeasible, got {model.PrimalRayLinObj}" \ No newline at end of file diff --git a/test/test_limit.py b/test/test_limit.py index 75ff6c1..0779380 100644 --- a/test/test_limit.py +++ b/test/test_limit.py @@ -12,28 +12,14 @@ # See the License for the specific language governing permissions and # limitations under the License. -# Copyright 2025 Haihao Lu -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. - -import time -import numpy as np -import scipy.sparse as sp -from cupdlpx import Model +import time +import numpy as np +import scipy.sparse as sp +from cupdlpx import Model, PDLP SEED = 42 -def test_time_limit(atol): +def test_time_limit(): """ Test time limit for large sparse LP. """ @@ -57,7 +43,7 @@ def test_time_limit(atol): tock = time.time() # check status assert hasattr(model, "Status"), "Model.Status not exposed." - assert model.Status == "TIME_LIMIT", f"Unexpected termination status: {model.Status}" + assert model.Status == PDLP.TIME_LIMIT, f"Unexpected termination status: {model.Status}" # check solving time solving_time = tock - tick assert solving_time < 1, f"Solving time exceeded limit a lot: {solving_time} seconds" @@ -65,7 +51,7 @@ def test_time_limit(atol): assert model.Runtime < 1, f"Internal solving time exceeded limit a lot: {model.Runtime} seconds" -def test_iters_limit(atol): +def test_iters_limit(): """ Test iteration limit for large sparse LP. """ @@ -80,14 +66,14 @@ def test_iters_limit(atol): ub = None model = Model(c, A, l, u, lb, ub) # turn off output - model.setParams(OutputFlag=False) + model.setParams(OutputFlag=False, Presolve=False) # set iteration limit model.setParams(IterationLimit=200) # optimize model.optimize() # check status assert hasattr(model, "Status"), "Model.Status not exposed." - assert model.Status == "ITERATION_LIMIT", f"Unexpected termination status: {model.Status}" + assert model.Status == PDLP.ITERATION_LIMIT, f"Unexpected termination status: {model.Status}" # check solving time assert hasattr(model, "IterCount"), "Model.IterCount not exposed." - assert model.IterCount <= 200, f"Internal iteration count exceeded limit a lot: {model.IterCount} seconds" \ No newline at end of file + assert model.IterCount <= 200, f"Internal iteration count exceeded limit a lot: {model.IterCount} seconds" diff --git a/test/test_matrix_formats.py b/test/test_matrix_formats.py index 74d97f7..f5e66db 100644 --- a/test/test_matrix_formats.py +++ b/test/test_matrix_formats.py @@ -13,7 +13,7 @@ # limitations under the License. import numpy as np -from cupdlpx import Model +from cupdlpx import Model, PDLP import scipy.sparse as sp def test_csr(base_lp_data, atol): @@ -30,7 +30,7 @@ def test_csr(base_lp_data, atol): model.optimize() # check status assert hasattr(model, "Status"), "Model.Status not exposed." - assert model.Status == "OPTIMAL", f"Unexpected termination status: {model.Status}" + assert model.Status == PDLP.OPTIMAL, f"Unexpected termination status: {model.Status}" # check primal solution assert hasattr(model, "X"), "Model.X (primal solution) not exposed." assert np.allclose(model.X, [1, 2], atol=atol), f"Unexpected primal solution: {model.X}" @@ -56,7 +56,7 @@ def test_csc(base_lp_data, atol): model.optimize() # check status assert hasattr(model, "Status"), "Model.Status not exposed." - assert model.Status == "OPTIMAL", f"Unexpected termination status: {model.Status}" + assert model.Status == PDLP.OPTIMAL, f"Unexpected termination status: {model.Status}" # check primal solution assert hasattr(model, "X"), "Model.X (primal solution) not exposed." assert np.allclose(model.X, [1, 2], atol=atol), f"Unexpected primal solution: {model.X}" @@ -82,7 +82,7 @@ def test_coo(base_lp_data, atol): model.optimize() # check status assert hasattr(model, "Status"), "Model.Status not exposed." - assert model.Status == "OPTIMAL", f"Unexpected termination status: {model.Status}" + assert model.Status == PDLP.OPTIMAL, f"Unexpected termination status: {model.Status}" # check primal solution assert hasattr(model, "X"), "Model.X (primal solution) not exposed." assert np.allclose(model.X, [1, 2], atol=atol), f"Unexpected primal solution: {model.X}" diff --git a/test/test_numerical.py b/test/test_numerical.py index 30415db..3f6d4f8 100644 --- a/test/test_numerical.py +++ b/test/test_numerical.py @@ -14,7 +14,7 @@ import numpy as np import scipy.sparse as sp -from cupdlpx import Model +from cupdlpx import Model, PDLP SEED = 42 @@ -38,7 +38,7 @@ def test_random_sparse_lp(atol): model.optimize() # check status assert hasattr(model, "Status"), "Model.Status not exposed." - assert model.Status == "OPTIMAL", f"Unexpected termination status: {model.Status}" + assert model.Status == PDLP.OPTIMAL, f"Unexpected termination status: {model.Status}" # primal objective obj_primal = c @ model.X # dual objective diff --git a/test/test_presolve.py b/test/test_presolve.py index d58c613..4a6d940 100644 --- a/test/test_presolve.py +++ b/test/test_presolve.py @@ -35,7 +35,7 @@ def test_presolve_as_optimal(base_lp_data, atol): model.optimize() # check status assert hasattr(model, "Status"), "Model.Status not exposed." - assert model.Status == "OPTIMAL", f"Unexpected termination status: {model.Status}" + assert model.Status == PDLP.OPTIMAL, f"Unexpected termination status: {model.Status}" # check primal solution assert hasattr(model, "X"), "Model.X (primal solution) not exposed." assert np.allclose(model.X, [1, 2], atol=atol), f"Unexpected primal solution: {model.X}" @@ -72,8 +72,7 @@ def test_presolve_as_infeasible(base_lp_data, atol): model.optimize() # check status assert hasattr(model, "Status"), "Model.Status not exposed." - assert model.Status == "PRIMAL_INFEASIBLE", f"Unexpected termination status: {model.Status}" - assert model.StatusCode == PDLP.PRIMAL_INFEASIBLE, f"Unexpected termination status code: {model.StatusCode}" + assert model.Status == PDLP.PRIMAL_INFEASIBLE, f"Unexpected termination status: {model.Status}" # check dual ray #assert model.DualRayObj > atol, f"DualRayObj should be positive for dual infeasible, got {model.DualRayObj}" diff --git a/test/test_read_mps.py b/test/test_read_mps.py index 91937c4..6167fa6 100644 --- a/test/test_read_mps.py +++ b/test/test_read_mps.py @@ -146,7 +146,7 @@ def test_read_and_optimize_minimize(mps_min_file, atol): model = cupdlpx.read(str(mps_min_file)) model.setParams(OutputFlag=False, Presolve=False) model.optimize() - assert model.Status == "OPTIMAL", f"Unexpected termination status: {model.Status}" + assert model.Status == PDLP.OPTIMAL, f"Unexpected termination status: {model.Status}" assert np.allclose(model.X, [1, 2], atol=atol), f"Unexpected primal solution: {model.X}" assert np.isclose(model.ObjVal, 3, atol=atol), f"Unexpected objective value: {model.ObjVal}" @@ -159,6 +159,6 @@ def test_read_and_optimize_maximize(mps_max_file, atol): model = cupdlpx.read(str(mps_max_file)) model.setParams(OutputFlag=False, Presolve=False) model.optimize() - assert model.Status == "OPTIMAL", f"Unexpected termination status: {model.Status}" + assert model.Status == PDLP.OPTIMAL, f"Unexpected termination status: {model.Status}" assert np.allclose(model.X, [1.5, 1.75], atol=atol), f"Unexpected primal solution: {model.X}" assert np.isclose(model.ObjVal, 3.25, atol=atol), f"Unexpected objective value: {model.ObjVal}" diff --git a/test/test_warm_start.py b/test/test_warm_start.py index af4335c..ea20565 100644 --- a/test/test_warm_start.py +++ b/test/test_warm_start.py @@ -13,7 +13,6 @@ # limitations under the License. import pytest -import warnings import numpy as np from cupdlpx import Model, PDLP @@ -39,7 +38,7 @@ def test_warm_start(base_lp_data, atol): # cold start baseline model.optimize() assert hasattr(model, "Status"), "Model.Status not exposed." - assert model.Status == "OPTIMAL", f"Unexpected termination status (cold): {model.Status}" + assert model.Status == PDLP.OPTIMAL, f"Unexpected termination status (cold): {model.Status}" assert hasattr(model, "IterCount"), "Model.IterCount not exposed." baseline_iters = model.IterCount # set warm start values @@ -48,7 +47,7 @@ def test_warm_start(base_lp_data, atol): model.optimize() # check status assert hasattr(model, "Status"), "Model.Status not exposed." - assert model.Status == "OPTIMAL", f"Unexpected termination status: {model.Status}" + assert model.Status == PDLP.OPTIMAL, f"Unexpected termination status: {model.Status}" # check primal solution assert hasattr(model, "X"), "Model.X (primal solution) not exposed." assert np.allclose(model.X, [1, 2], atol=atol), f"Unexpected primal solution: {model.X}" @@ -89,7 +88,7 @@ def test_warm_start_primal(base_lp_data, atol): model.optimize() # check status assert hasattr(model, "Status"), "Model.Status not exposed." - assert model.Status == "OPTIMAL", f"Unexpected termination status: {model.Status}" + assert model.Status == PDLP.OPTIMAL, f"Unexpected termination status: {model.Status}" # check primal solution assert hasattr(model, "X"), "Model.X (primal solution) not exposed." assert np.allclose(model.X, [1, 2], atol=atol), f"Unexpected primal solution: {model.X}" @@ -125,7 +124,7 @@ def test_warm_start_dual(base_lp_data, atol): model.optimize() # check status assert hasattr(model, "Status"), "Model.Status not exposed." - assert model.Status == "OPTIMAL", f"Unexpected termination status: {model.Status}" + assert model.Status == PDLP.OPTIMAL, f"Unexpected termination status: {model.Status}" # check primal solution assert hasattr(model, "X"), "Model.X (primal solution) not exposed." assert np.allclose(model.X, [1, 2], atol=atol), f"Unexpected primal solution: {model.X}" @@ -156,7 +155,7 @@ def test_clear_warm_start(base_lp_data, atol): model.optimize() # check status assert hasattr(model, "Status"), "Model.Status not exposed." - assert model.Status == "OPTIMAL", f"Unexpected termination status: {model.Status}" + assert model.Status == PDLP.OPTIMAL, f"Unexpected termination status: {model.Status}" # check primal solution assert hasattr(model, "X"), "Model.X (primal solution) not exposed." assert np.allclose(model.X, [1, 2], atol=atol), f"Unexpected primal solution: {model.X}" @@ -167,29 +166,33 @@ def test_clear_warm_start(base_lp_data, atol): assert hasattr(model, "ObjVal"), "Model.ObjVal (objective value) not exposed." assert np.isclose(model.ObjVal, 3, atol=atol), f"Unexpected objective value: {model.ObjVal}" -def test_warm_start_wrong_size_fallback(base_lp_data, atol): +def test_warm_start_wrong_size_raises(base_lp_data): """ - Verify that warm start with wrong size falls back to cold start with a warning. + Verify that warm start values with wrong sizes are rejected with ValueError. """ # setup model c, A, l, u, lb, ub = base_lp_data model = Model(c, A, l, u, lb, ub) # turn off output model.setParams(OutputFlag=False, Presolve=False) - # set warm start values with wrong size - with pytest.warns(RuntimeWarning): - model.setWarmStart(primal=[1], dual=[1, 1]) # wrong sizes - # optimize - model.optimize() - # check status - assert hasattr(model, "Status"), "Model.Status not exposed." - assert model.Status == "OPTIMAL", f"Unexpected termination status: {model.Status}" - # check primal solution - assert hasattr(model, "X"), "Model.X (primal solution) not exposed." - assert np.allclose(model.X, [1, 2], atol=atol), f"Unexpected primal solution: {model.X}" - # check dual solution - assert hasattr(model, "Pi"), "Model.Pi (dual solution) not exposed." - assert np.allclose(model.Pi, [1, -1, 0], atol=atol), f"Unexpected dual solution: {model.Pi}" - # check objective - assert hasattr(model, "ObjVal"), "Model.ObjVal (objective value) not exposed." - assert np.isclose(model.ObjVal, 3, atol=atol), f"Unexpected objective value: {model.ObjVal}" \ No newline at end of file + # wrong-size warm start values are rejected + with pytest.raises(ValueError): + model.setWarmStart(primal=[1]) # wrong size + with pytest.raises(ValueError): + model.setWarmStart(dual=[1, 1]) # wrong size + # rejected values must not be stored + assert model._primal_start is None, "Rejected primal warm start was stored." + assert model._dual_start is None, "Rejected dual warm start was stored." + + +def test_warm_start_update_is_transactional(base_lp_data): + c, A, l, u, lb, ub = base_lp_data + model = Model(c, A, l, u, lb, ub) + model.setParams(OutputFlag=False, Presolve=False) + model.setWarmStart(primal=[1, 2], dual=[1, -1, 0]) + + with pytest.raises(ValueError): + model.setWarmStart(primal=[2, 1], dual=[1, 2]) + + assert np.allclose(model._primal_start, [1, 2]) + assert np.allclose(model._dual_start, [1, -1, 0])