From 22021180f63738082dd9d5b2be14f3d42c141519 Mon Sep 17 00:00:00 2001 From: mal84 Date: Wed, 29 Oct 2025 15:48:03 +0000 Subject: [PATCH 1/6] Implement interface to cupdlpx solver and required IO --- linopy/io.py | 46 +++++++++ linopy/model.py | 3 + linopy/solvers.py | 255 ++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 304 insertions(+) diff --git a/linopy/io.py b/linopy/io.py index 7065adbb..f74d83b7 100644 --- a/linopy/io.py +++ b/linopy/io.py @@ -27,6 +27,7 @@ from linopy.objective import Objective if TYPE_CHECKING: + from cupdlpx import Model as cupdlpxModel from highspy.highs import Highs from linopy.model import Model @@ -756,6 +757,51 @@ def to_highspy(m: Model, explicit_coordinate_names: bool = False) -> Highs: return h +def to_cupdlpx(m: Model, explicit_coordinate_names: bool = False) -> cupdlpxModel: + """ + Export the model to cupdlpx. + + This function does not write the model to intermediate files but directly + passes it to cupdlpx. + + cuPDLPx does not support named variables and constraints, so names + are not tracked by this function. + + Parameters + ---------- + m : linopy.Model + + Returns + ------- + model : cupdlpx.Model + """ + import cupdlpx + + # build model using canonical form matrices and vectors + # see https://github.com/MIT-Lu-Lab/cuPDLPx/tree/main/python#modeling + M = m.matrices + A = M.A.tocsr() # cuDPLPx only support CSR sparse matrix format + # linopy stores constraints as Ax ?= b and keeps track of inequality + # sense in M.sense. Convert to separate lower and upper bound vectors. + l = np.where(M.sense == ">", M.b, -np.inf) + u = np.where(M.sense == "<", M.b, np.inf) + + cu_model = cupdlpx.Model( + objective_vector=M.c, + constraint_matrix=A, + constraint_lower_bound=l, + constraint_upper_bound=u, + variable_lower_bound=M.lb, + variable_upper_bound=M.ub, + ) + + # change objective sense + if m.objective.sense == "max": + cu_model.ModelSense = cupdlpx.PDLP.MAXIMIZE + + return cu_model + + def to_block_files(m: Model, fn: Path) -> None: """ Write out the linopy model to a block structured output. diff --git a/linopy/model.py b/linopy/model.py index 3982b84d..79fbf2b4 100644 --- a/linopy/model.py +++ b/linopy/model.py @@ -50,6 +50,7 @@ ) from linopy.io import ( to_block_files, + to_cupdlpx, to_file, to_gurobipy, to_highspy, @@ -1510,4 +1511,6 @@ def reset_solution(self) -> None: to_highspy = to_highspy + to_cupdlpx = to_cupdlpx + to_block_files = to_block_files diff --git a/linopy/solvers.py b/linopy/solvers.py index 087ec75e..675a60c2 100644 --- a/linopy/solvers.py +++ b/linopy/solvers.py @@ -24,6 +24,7 @@ import pandas as pd from packaging.version import parse as parse_version +import linopy.io from linopy.constants import ( Result, Solution, @@ -59,6 +60,7 @@ "scip", "copt", "mindopt", + "cupdlpx", ] FILE_IO_APIS = ["lp", "lp-polars", "mps"] @@ -134,6 +136,16 @@ except coptpy.CoptError: pass +with contextlib.suppress(ModuleNotFoundError): + import cupdlpx + + try: + cupdlpx.Model(np.array([0.0]), np.array([[0.0]]), None, None) + available_solvers.append("cupdlpx") + except ImportError: + pass + + quadratic_solvers = [s for s in QUADRATIC_SOLVERS if s in available_solvers] logger = logging.getLogger(__name__) @@ -165,6 +177,7 @@ class SolverName(enum.Enum): COPT = "copt" MindOpt = "mindopt" PIPS = "pips" + cuPDLPx = "cupdlpx" def path_to_string(path: Path) -> str: @@ -2261,3 +2274,245 @@ def __init__( super().__init__(**solver_options) msg = "The PIPS solver interface is not yet implemented." raise NotImplementedError(msg) + + +class cuPDLPx(Solver[None]): + """ + Solver subclass for the cuPDLPx solver. cuPDLPx must be installed + with working GPU support for usage. Find the installation instructions + at https://github.com/MIT-Lu-Lab/cuPDLPx. + + The full list of solver options provided with the python interface + is documented at https://github.com/MIT-Lu-Lab/cuPDLPx/tree/main/python. + + Some example options are: + * LogToConsole : False by default. + * TimeLimit : 3600.0 by default. + * IterationLimit : 2147483647 by default. + + Attributes + ---------- + **solver_options + options for the given solver + """ + + def __init__( + self, + **solver_options: Any, + ) -> None: + super().__init__(**solver_options) + + def solve_problem_from_file( + self, + problem_fn: Path, + solution_fn: Path | None = None, + log_fn: Path | None = None, + warmstart_fn: Path | None = None, + basis_fn: Path | None = None, + env: EnvType | None = None, + ) -> Result: + """ + Solve a linear problem from a problem file using the solver cuPDLPx. + cuPDLPx does not currently support its own file IO, so this function + reads the problem file using linopy (only support netcf files) and + then passes the model to cuPDLPx for solving. + If the solution is feasible the function returns the + objective, solution and dual constraint variables. + + Parameters + ---------- + problem_fn : Path + Path to the problem file. + solution_fn : Path, optional + Path to the solution file. + log_fn : Path, optional + Path to the log file. + warmstart_fn : Path, optional + Path to the warmstart file. + basis_fn : Path, optional + Path to the basis file. + env : None, optional + Environment for the solver + + Returns + ------- + Result + """ + logger.warning( + "cuPDLPx doesn't currently support file IO. Building model from file using linopy." + ) + problem_fn_ = path_to_string(problem_fn) + + if problem_fn_.endswith(".netcdf"): + model: Model = linopy.io.read_netcdf(problem_fn_) + else: + msg = "linopy currently only supports reading models from netcdf files. Try using io_api='direct' instead." + raise NotImplementedError(msg) + + return self.solve_problem_from_model( + model, + solution_fn=solution_fn, + log_fn=log_fn, + warmstart_fn=warmstart_fn, + basis_fn=basis_fn, + env=env, + ) + + def solve_problem_from_model( + self, + model: Model, + solution_fn: Path | None = None, + log_fn: Path | None = None, + warmstart_fn: Path | None = None, + basis_fn: Path | None = None, + env: EnvType | None = None, + explicit_coordinate_names: bool = False, + ) -> Result: + """ + Solve a linear problem directly from a linopy model using the solver cuPDLPx. + If the solution is feasible the function returns the + objective, solution and dual constraint variables. + + Parameters + ---------- + model : linopy.model + Linopy model for the problem. + solution_fn : Path, optional + Path to the solution file. + log_fn : Path, optional + Path to the log file. + warmstart_fn : Path, optional + Path to the warmstart file. + basis_fn : Path, optional + Path to the basis file. + env : None, optional + Environment for the solver + explicit_coordinate_names : bool, optional + Transfer variable and constraint names to the solver (default: False) + + Returns + ------- + Result + """ + + if model.type in ["QP", "MILP"]: + msg = "cuPDLPx does not currently support QP or MILP problems." + raise NotImplementedError(msg) + + cu_model = model.to_cupdlpx() + + return self._solve( + cu_model, + l_model=model, + solution_fn=solution_fn, + log_fn=log_fn, + warmstart_fn=warmstart_fn, + basis_fn=basis_fn, + io_api="direct", + sense=model.sense, + ) + + def _solve( + self, + cu_model: cupdlpx.Model, + l_model: Model | None = None, + solution_fn: Path | None = None, + log_fn: Path | None = None, + warmstart_fn: Path | None = None, + basis_fn: Path | None = None, + io_api: str | None = None, + sense: str | None = None, + ) -> Result: + """ + Solve a linear problem from a cupdlpx.Model object. + + Parameters + ---------- + cu_model: cupdlpx.Model + cupdlpx object. + solution_fn : Path, optional + Path to the solution file. + log_fn : Path, optional + Path to the log file. + warmstart_fn : Path, optional + Path to the warmstart file. + basis_fn : Path, optional + Path to the basis file. + model : linopy.model, optional + Linopy model for the problem. + io_api: str + io_api of the problem. For direct API from linopy model this is "direct". + sense: str + "min" or "max" + + Returns + ------- + Result + """ + + # see https://github.com/MIT-Lu-Lab/cuPDLPx/blob/main/python/cupdlpx/PDLP.py + CONDITION_MAP: dict[int, TerminationCondition] = { + cupdlpx.PDLP.OPTIMAL: TerminationCondition.optimal, + cupdlpx.PDLP.PRIMAL_INFEASIBLE: TerminationCondition.infeasible, + cupdlpx.PDLP.DUAL_INFEASIBLE: TerminationCondition.infeasible_or_unbounded, + cupdlpx.PDLP.TIME_LIMIT: TerminationCondition.time_limit, + cupdlpx.PDLP.ITERATION_LIMIT: TerminationCondition.iteration_limit, + cupdlpx.PDLP.UNSPECIFIED: TerminationCondition.unknown, + } + + self._set_solver_params(cu_model) + + if warmstart_fn is not None: + # cuPDLPx supports warmstart, but there currently isn't the tooling + # to read it in from a file + raise NotImplementedError("Warmstarting not yet implemented for cuPDLPx.") + else: + cu_model.clearWarmStart() + + if basis_fn is not None: + logger.warning("Basis files are not supported by cuPDLPx. Ignoring.") + + # solve + cu_model.optimize() + + # parse solution and output + if solution_fn is not None: + raise NotImplementedError( + "Solution file output not yet implemented for cuPDLPx." + ) + + termination_condition = CONDITION_MAP.get( + cu_model.StatusCode, cu_model.StatusCode + ) + status = Status.from_termination_condition(termination_condition) + status.legacy_status = cu_model.Status # cuPDLPx status message + + def get_solver_solution() -> Solution: + objective = cu_model.ObjVal + + vlabels = None if l_model is None else l_model.matrices.vlabels + clabels = None if l_model is None else l_model.matrices.clabels + + sol = pd.Series(cu_model.X, vlabels, dtype=float) + dual = pd.Series(cu_model.Pi, clabels, dtype=float) + + if cu_model.ModelSense == cupdlpx.PDLP.MAXIMIZE: + dual *= -1 # flip sign of duals for max problems + + return Solution(sol, dual, objective) + + solution = self.safe_get_solution(status=status, func=get_solver_solution) + solution = maybe_adjust_objective_sign(solution, io_api, sense) + + # see https://github.com/MIT-Lu-Lab/cuPDLPx/tree/main/python#solution-attributes + return Result(status, solution, cu_model) + + def _set_solver_params(self, cu_model: cupdlpx.Model) -> None: + """ + Set solver options for cuPDLPx model. + + For list of available options, see + https://github.com/MIT-Lu-Lab/cuPDLPx/tree/main/python#parameters + """ + for k, v in self.solver_options.items(): + cu_model.setParam(k, v) From 7b69035b62ae5850e6cb55efd2df5d500c5b9dd3 Mon Sep 17 00:00:00 2001 From: mal84 Date: Wed, 29 Oct 2025 15:48:22 +0000 Subject: [PATCH 2/6] Update tests --- test/test_optimization.py | 204 ++++++++++++++++++++++---------------- test/test_solvers.py | 3 + 2 files changed, 123 insertions(+), 84 deletions(-) diff --git a/test/test_optimization.py b/test/test_optimization.py index da1050fc..b8bcd8c9 100644 --- a/test/test_optimization.py +++ b/test/test_optimization.py @@ -16,7 +16,7 @@ import pandas as pd import pytest import xarray as xr -from xarray.testing import assert_equal +from xarray.testing import assert_allclose, assert_equal from linopy import GREATER_EQUAL, LESS_EQUAL, Model, solvers from linopy.common import to_path @@ -32,12 +32,12 @@ # mps io is only supported via highspy io_apis.append("mps") - +file_io_solvers = [s for s in available_solvers if s not in ["cupdlpx"]] params: list[tuple[str, str, bool]] = list( - itertools.product(available_solvers, io_apis, explicit_coordinate_names) + itertools.product(file_io_solvers, io_apis, explicit_coordinate_names) ) -direct_solvers: list[str] = ["gurobi", "highs", "mosek"] +direct_solvers: list[str] = ["gurobi", "highs", "mosek", "cupdlpx"] for solver in direct_solvers: if solver in available_solvers: params.append((solver, "direct", False)) @@ -53,6 +53,9 @@ if platform.system() == "Windows" and "scip" in feasible_quadratic_solvers: feasible_quadratic_solvers.remove("scip") +feasible_mip_solvers: list[str] = available_solvers.copy() +feasible_mip_solvers.remove("cupdlpx") # cuPDLPx does not support MIP yet + def test_print_solvers(capsys: Any) -> None: with capsys.disabled(): @@ -407,7 +410,7 @@ def test_anonymous_constraint( model.solve( solver, io_api=io_api, explicit_coordinate_names=explicit_coordinate_names ) - assert_equal(model.solution, model_anonymous_constraint.solution) + assert_allclose(model.solution, model_anonymous_constraint.solution) @pytest.mark.parametrize("solver,io_api,explicit_coordinate_names", params) @@ -475,6 +478,7 @@ def test_solver_time_limit_options( "mosek": {"MSK_DPAR_OPTIMIZER_MAX_TIME": 1}, "mindopt": {"MaxTime": 1}, "copt": {"TimeLimit": 1}, + "cupdlpx": {"TimeLimit": 1}, } status, condition = model.solve( solver, @@ -508,7 +512,7 @@ def test_duplicated_variables( ) -> None: status, condition = model_with_duplicated_variables.solve(solver, io_api=io_api) assert status == "ok" - assert all(model_with_duplicated_variables.solution["x"] == 5) + assert all(np.isclose(model_with_duplicated_variables.solution["x"], 5, rtol=1e-4)) @pytest.mark.parametrize("solver,io_api,explicit_coordinate_names", params) @@ -523,9 +527,15 @@ def test_non_aligned_variables( ) assert status == "ok" with pytest.warns(UserWarning): - assert model_with_non_aligned_variables.solution["x"][0] == 0 - assert model_with_non_aligned_variables.solution["x"][-1] == 10.5 - assert model_with_non_aligned_variables.solution["y"][0] == 10.5 + assert np.isclose( + model_with_non_aligned_variables.solution["x"][0], 0, rtol=1e-4 + ) + assert np.isclose( + model_with_non_aligned_variables.solution["x"][-1], 10.5, rtol=1e-4 + ) + assert np.isclose( + model_with_non_aligned_variables.solution["y"][0], 10.5, rtol=1e-4 + ) assert np.isnan(model_with_non_aligned_variables.solution["y"][-1]) for dtype in model_with_non_aligned_variables.solution.dtypes.values(): @@ -540,16 +550,17 @@ def test_set_files( io_api: str, explicit_coordinate_names: bool, ) -> None: - status, condition = model.solve( - solver, - io_api=io_api, - explicit_coordinate_names=explicit_coordinate_names, - problem_fn=tmp_path / "problem.lp", - solution_fn=tmp_path / "solution.sol", - log_fn=tmp_path / "logging.log", - keep_files=False, - ) - assert status == "ok" + if solver in file_io_solvers: + status, condition = model.solve( + solver, + io_api=io_api, + explicit_coordinate_names=explicit_coordinate_names, + problem_fn=tmp_path / "problem.lp", + solution_fn=tmp_path / "solution.sol", + log_fn=tmp_path / "logging.log", + keep_files=False, + ) + assert status == "ok" @pytest.mark.parametrize("solver,io_api,explicit_coordinate_names", params) @@ -560,26 +571,33 @@ def test_set_files_and_keep_files( io_api: str, explicit_coordinate_names: bool, ) -> None: - status, condition = model.solve( - solver, - io_api=io_api, - explicit_coordinate_names=explicit_coordinate_names, - problem_fn=tmp_path / "problem.lp", - solution_fn=tmp_path / "solution.sol", - log_fn=tmp_path / "logging.log", - keep_files=True, - ) - assert status == "ok" - if io_api != "direct" and solver != "xpress": - assert (tmp_path / "problem.lp").exists() - assert (tmp_path / "solution.sol").exists() - assert (tmp_path / "logging.log").exists() + if solver in file_io_solvers: + status, condition = model.solve( + solver, + io_api=io_api, + explicit_coordinate_names=explicit_coordinate_names, + problem_fn=tmp_path / "problem.lp", + solution_fn=tmp_path / "solution.sol", + log_fn=tmp_path / "logging.log", + keep_files=True, + ) + assert status == "ok" + if io_api != "direct" and solver != "xpress": + assert (tmp_path / "problem.lp").exists() + assert (tmp_path / "solution.sol").exists() + assert (tmp_path / "logging.log").exists() @pytest.mark.parametrize("solver,io_api,explicit_coordinate_names", params) def test_infeasible_model( model: Model, solver: str, io_api: str, explicit_coordinate_names: bool ) -> None: + if solver == "cupdlpx": + pytest.skip( + "Ongoing issue with cuPDLPx causes it to hang for some unbounded problems. " + "See https://github.com/MIT-Lu-Lab/cuPDLPx/issues/9." + ) + model.add_constraints([(1, "x")], "<=", 0) model.add_constraints([(1, "y")], "<=", 0) @@ -604,10 +622,11 @@ def test_infeasible_model( def test_model_with_inf( model_with_inf: Model, solver: str, io_api: str, explicit_coordinate_names: bool ) -> None: - status, condition = model_with_inf.solve(solver, io_api=io_api) - assert condition == "optimal" - assert (model_with_inf.solution.x == 0).all() - assert (model_with_inf.solution.y == 10).all() + if solver in feasible_mip_solvers: + status, condition = model_with_inf.solve(solver, io_api=io_api) + assert condition == "optimal" + assert (model_with_inf.solution.x == 0).all() + assert (model_with_inf.solution.y == 10).all() @pytest.mark.parametrize( @@ -617,13 +636,14 @@ def test_model_with_inf( def test_milp_binary_model( milp_binary_model: Model, solver: str, io_api: str, explicit_coordinate_names: bool ) -> None: - status, condition = milp_binary_model.solve( - solver, io_api=io_api, explicit_coordinate_names=explicit_coordinate_names - ) - assert condition == "optimal" - assert ( - (milp_binary_model.solution.y == 1) | (milp_binary_model.solution.y == 0) - ).all() + if solver in feasible_mip_solvers: + status, condition = milp_binary_model.solve( + solver, io_api=io_api, explicit_coordinate_names=explicit_coordinate_names + ) + assert condition == "optimal" + assert ( + (milp_binary_model.solution.y == 1) | (milp_binary_model.solution.y == 0) + ).all() @pytest.mark.parametrize( @@ -636,13 +656,15 @@ def test_milp_binary_model_r( io_api: str, explicit_coordinate_names: bool, ) -> None: - status, condition = milp_binary_model_r.solve( - solver, io_api=io_api, explicit_coordinate_names=explicit_coordinate_names - ) - assert condition == "optimal" - assert ( - (milp_binary_model_r.solution.x == 1) | (milp_binary_model_r.solution.x == 0) - ).all() + if solver in feasible_mip_solvers: + status, condition = milp_binary_model_r.solve( + solver, io_api=io_api, explicit_coordinate_names=explicit_coordinate_names + ) + assert condition == "optimal" + assert ( + (milp_binary_model_r.solution.x == 1) + | (milp_binary_model_r.solution.x == 0) + ).all() @pytest.mark.parametrize( @@ -652,11 +674,12 @@ def test_milp_binary_model_r( def test_milp_model( milp_model: Model, solver: str, io_api: str, explicit_coordinate_names: bool ) -> None: - status, condition = milp_model.solve( - solver, io_api=io_api, explicit_coordinate_names=explicit_coordinate_names - ) - assert condition == "optimal" - assert ((milp_model.solution.y == 9) | (milp_model.solution.x == 0.5)).all() + if solver in feasible_mip_solvers: + status, condition = milp_model.solve( + solver, io_api=io_api, explicit_coordinate_names=explicit_coordinate_names + ) + assert condition == "optimal" + assert ((milp_model.solution.y == 9) | (milp_model.solution.x == 0.5)).all() @pytest.mark.parametrize( @@ -666,14 +689,19 @@ def test_milp_model( def test_milp_model_r( milp_model_r: Model, solver: str, io_api: str, explicit_coordinate_names: bool ) -> None: - # MPS format by Highs wrong, see https://github.com/ERGO-Code/HiGHS/issues/1325 - # skip it - if io_api != "mps": - status, condition = milp_model_r.solve( - solver, io_api=io_api, explicit_coordinate_names=explicit_coordinate_names - ) - assert condition == "optimal" - assert ((milp_model_r.solution.x == 11) | (milp_model_r.solution.y == 0)).all() + if solver in feasible_mip_solvers: + # MPS format by Highs wrong, see https://github.com/ERGO-Code/HiGHS/issues/1325 + # skip it + if io_api != "mps": + status, condition = milp_model_r.solve( + solver, + io_api=io_api, + explicit_coordinate_names=explicit_coordinate_names, + ) + assert condition == "optimal" + assert ( + (milp_model_r.solution.x == 11) | (milp_model_r.solution.y == 0) + ).all() @pytest.mark.parametrize( @@ -798,13 +826,13 @@ def test_quadratic_model_unbounded( def test_modified_model( modified_model: Model, solver: str, io_api: str, explicit_coordinate_names: bool ) -> None: - status, condition = modified_model.solve( - solver, io_api=io_api, explicit_coordinate_names=explicit_coordinate_names - ) - - assert condition == "optimal" - assert (modified_model.solution.x == 0).all() - assert (modified_model.solution.y == 10).all() + if solver in feasible_mip_solvers: + status, condition = modified_model.solve( + solver, io_api=io_api, explicit_coordinate_names=explicit_coordinate_names + ) + assert condition == "optimal" + assert (modified_model.solution.x == 0).all() + assert (modified_model.solution.y == 10).all() @pytest.mark.parametrize("solver,io_api,explicit_coordinate_names", params) @@ -822,7 +850,7 @@ def test_masked_variable_model( assert y.solution[-2:].isnull().all() assert y.solution[:-2].notnull().all() assert x.solution.notnull().all() - assert (x.solution[-2:] == 10).all() + assert (np.isclose(x.solution[-2:], 10, rtol=2.5e-4)).all() # Squeeze in solution getter for expressions with masked variables assert_equal(x.add(y).solution, x.solution + y.solution.fillna(0)) @@ -837,8 +865,8 @@ def test_masked_constraint_model( masked_constraint_model.solve( solver, io_api=io_api, explicit_coordinate_names=explicit_coordinate_names ) - assert (masked_constraint_model.solution.y[:-2] == 10).all() - assert (masked_constraint_model.solution.y[-2:] == 5).all() + assert (np.isclose(masked_constraint_model.solution.y[:-2], 10, rtol=1e-4)).all() + assert (np.isclose(masked_constraint_model.solution.y[-2:], 5, rtol=2e-4)).all() @pytest.mark.parametrize("solver,io_api,explicit_coordinate_names", params) @@ -849,6 +877,9 @@ def test_basis_and_warmstart( io_api: str, explicit_coordinate_names: bool, ) -> None: + if solver == "cupdlpx": + pytest.skip("cuPDLPx does not yet support warmstart in the Python API.") + basis_fn = tmp_path / "basis.bas" model.solve( solver, @@ -872,14 +903,15 @@ def test_solution_fn_parent_dir_doesnt_exist( explicit_coordinate_names: bool, tmp_path: Any, ) -> None: - solution_fn = tmp_path / "non_existent_dir" / "non_existent_file" - status, condition = model.solve( - solver, - io_api=io_api, - explicit_coordinate_names=explicit_coordinate_names, - solution_fn=solution_fn, - ) - assert status == "ok" + if solver in file_io_solvers: + solution_fn = tmp_path / "non_existent_dir" / "non_existent_file" + status, condition = model.solve( + solver, + io_api=io_api, + explicit_coordinate_names=explicit_coordinate_names, + solution_fn=solution_fn, + ) + assert status == "ok" @pytest.mark.parametrize("solver", available_solvers) @@ -892,7 +924,7 @@ def test_non_supported_solver_io(model: Model, solver: str) -> None: def test_solver_attribute_getter( model: Model, solver: str, io_api: str, explicit_coordinate_names: bool ) -> None: - model.solve(solver) + model.solve(solver, io_api=io_api) if solver != "gurobi": with pytest.raises(NotImplementedError): model.variables.get_solver_attribute("RC") @@ -921,7 +953,11 @@ def test_model_resolve( ) assert status == "ok" # x = -0.75, y = 3.0 - assert np.isclose(model.objective.value or 0, 5.25) + + if solver == "cupdlpx": # this solver has low resolution + assert np.isclose(model.objective.value or 0, 5.25, rtol=1e-3) + else: + assert np.isclose(model.objective.value or 0, 5.25) @pytest.mark.parametrize( diff --git a/test/test_solvers.py b/test/test_solvers.py index 129c1e0b..09ff5651 100644 --- a/test/test_solvers.py +++ b/test/test_solvers.py @@ -53,6 +53,9 @@ def test_free_mps_solution_parsing(solver: str, tmp_path: Path) -> None: except ValueError: raise ValueError(f"Solver '{solver}' is not recognized") + if solver_class == solvers.cuPDLPx: + pytest.skip("cuPDLPx does not currently support file IO.") + # Write the MPS file to the temporary directory mps_file = tmp_path / "problem.mps" mps_file.write_text(free_mps_problem) From b683654e10934ba8be7b3cc3d75a9e3edd2968c3 Mon Sep 17 00:00:00 2001 From: mal84 <89249742+mal84emma@users.noreply.github.com> Date: Wed, 29 Oct 2025 15:56:22 +0000 Subject: [PATCH 3/6] Add to docs & add package dependency --- README.md | 1 + pyproject.toml | 7 +++---- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/README.md b/README.md index 3cd19df8..644b556c 100644 --- a/README.md +++ b/README.md @@ -149,6 +149,7 @@ Fri 0 4 * [Cplex](https://www.ibm.com/de-de/analytics/cplex-optimizer) * [MOSEK](https://www.mosek.com/) * [COPT](https://www.shanshu.ai/copt) +* [cuPDLPx](https://github.com/MIT-Lu-Lab/cuPDLPx) Note that these do have to be installed by the user separately. diff --git a/pyproject.toml b/pyproject.toml index b5105230..b4922ca9 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -37,7 +37,7 @@ dependencies = [ "tqdm", "deprecation", "google-cloud-storage", - "requests" + "requests", ] [project.urls] @@ -81,6 +81,7 @@ solvers = [ "coptpy!=7.2.1", "xpress; platform_system != 'Darwin' and python_version < '3.11'", "pyscipopt; platform_system != 'Darwin'", + "cupdlpx>=0.1.2", ] [tool.setuptools.packages.find] @@ -102,9 +103,7 @@ branch = true source = ["linopy"] omit = ["test/*"] [tool.coverage.report] -exclude_also = [ - "if TYPE_CHECKING:", -] +exclude_also = ["if TYPE_CHECKING:"] [tool.mypy] exclude = ['dev/*', 'examples/*', 'benchmark/*', 'doc/*'] From 6416b7189c66f77bb5b6d943e8309a0565265586 Mon Sep 17 00:00:00 2001 From: mal84 <89249742+mal84emma@users.noreply.github.com> Date: Wed, 29 Oct 2025 16:02:49 +0000 Subject: [PATCH 4/6] Add comment for release notes --- doc/release_notes.rst | 1 + 1 file changed, 1 insertion(+) diff --git a/doc/release_notes.rst b/doc/release_notes.rst index 079330eb..e4bd8c92 100644 --- a/doc/release_notes.rst +++ b/doc/release_notes.rst @@ -11,6 +11,7 @@ Upcoming Version * Harmonize dtypes before concatenation in lp file writing to avoid dtype mismatch errors. This error occurred when creating and storing models in netcdf format using windows machines and loading and solving them on linux machines. * Add option to use polars series as constant input * Fix expression merge to explicitly use outer join when combining expressions with disjoint coordinates for consistent behavior across xarray versions +* Add support for GPU-accelerated solver [cuPDLPx](https://github.com/MIT-Lu-Lab/cuPDLPx) Version 0.5.6 -------------- From e22f55ca6b03e5a86f7e7adb08da5cb3691ffd39 Mon Sep 17 00:00:00 2001 From: mal84 <89249742+mal84emma@users.noreply.github.com> Date: Thu, 30 Oct 2025 11:59:37 +0000 Subject: [PATCH 5/6] Correction: add support for equality constraints (same l & u bounds) --- linopy/io.py | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/linopy/io.py b/linopy/io.py index f74d83b7..6d23c5af 100644 --- a/linopy/io.py +++ b/linopy/io.py @@ -783,8 +783,16 @@ def to_cupdlpx(m: Model, explicit_coordinate_names: bool = False) -> cupdlpxMode A = M.A.tocsr() # cuDPLPx only support CSR sparse matrix format # linopy stores constraints as Ax ?= b and keeps track of inequality # sense in M.sense. Convert to separate lower and upper bound vectors. - l = np.where(M.sense == ">", M.b, -np.inf) - u = np.where(M.sense == "<", M.b, np.inf) + l = np.where( + np.logical_or(np.equal(M.sense, ">"), np.equal(M.sense, "=")), + M.b, + -np.inf, + ) + u = np.where( + np.logical_or(np.equal(M.sense, "<"), np.equal(M.sense, "=")), + M.b, + np.inf, + ) cu_model = cupdlpx.Model( objective_vector=M.c, From 81b202be0eb4355006b4aa8eab44155a980ea7f3 Mon Sep 17 00:00:00 2001 From: mal84 <89249742+mal84emma@users.noreply.github.com> Date: Fri, 31 Oct 2025 16:07:33 +0000 Subject: [PATCH 6/6] Parameterize solution tolerance in tests to allow different standard for CPU and GPU solver precisions --- test/test_optimization.py | 32 ++++++++++++++++++++------------ 1 file changed, 20 insertions(+), 12 deletions(-) diff --git a/test/test_optimization.py b/test/test_optimization.py index b8bcd8c9..0a0ea156 100644 --- a/test/test_optimization.py +++ b/test/test_optimization.py @@ -46,7 +46,6 @@ params.append(("mosek", "lp", False)) params.append(("mosek", "lp", True)) - feasible_quadratic_solvers: list[str] = quadratic_solvers # There seems to be a bug in scipopt with quadratic models on windows, see # https://github.com/PyPSA/linopy/actions/runs/7615240686/job/20739454099?pr=78 @@ -56,6 +55,11 @@ feasible_mip_solvers: list[str] = available_solvers.copy() feasible_mip_solvers.remove("cupdlpx") # cuPDLPx does not support MIP yet +gpu_solvers: list[str] = ["cupdlpx"] + +CPU_SOL_TOL: float = 1e-5 # numpy default +GPU_SOL_TOL: float = 2.5e-4 # gpu solvers typically have lower numerical precision + def test_print_solvers(capsys: Any) -> None: with capsys.disabled(): @@ -512,7 +516,8 @@ def test_duplicated_variables( ) -> None: status, condition = model_with_duplicated_variables.solve(solver, io_api=io_api) assert status == "ok" - assert all(np.isclose(model_with_duplicated_variables.solution["x"], 5, rtol=1e-4)) + tol = GPU_SOL_TOL if solver in gpu_solvers else CPU_SOL_TOL + assert all(np.isclose(model_with_duplicated_variables.solution["x"], 5, rtol=tol)) @pytest.mark.parametrize("solver,io_api,explicit_coordinate_names", params) @@ -526,15 +531,18 @@ def test_non_aligned_variables( solver, io_api=io_api, explicit_coordinate_names=explicit_coordinate_names ) assert status == "ok" + + tol = GPU_SOL_TOL if solver in gpu_solvers else CPU_SOL_TOL + with pytest.warns(UserWarning): assert np.isclose( - model_with_non_aligned_variables.solution["x"][0], 0, rtol=1e-4 + model_with_non_aligned_variables.solution["x"][0], 0, rtol=tol ) assert np.isclose( - model_with_non_aligned_variables.solution["x"][-1], 10.5, rtol=1e-4 + model_with_non_aligned_variables.solution["x"][-1], 10.5, rtol=tol ) assert np.isclose( - model_with_non_aligned_variables.solution["y"][0], 10.5, rtol=1e-4 + model_with_non_aligned_variables.solution["y"][0], 10.5, rtol=tol ) assert np.isnan(model_with_non_aligned_variables.solution["y"][-1]) @@ -850,7 +858,8 @@ def test_masked_variable_model( assert y.solution[-2:].isnull().all() assert y.solution[:-2].notnull().all() assert x.solution.notnull().all() - assert (np.isclose(x.solution[-2:], 10, rtol=2.5e-4)).all() + tol = GPU_SOL_TOL if solver in gpu_solvers else CPU_SOL_TOL + assert (np.isclose(x.solution[-2:], 10, rtol=tol)).all() # Squeeze in solution getter for expressions with masked variables assert_equal(x.add(y).solution, x.solution + y.solution.fillna(0)) @@ -865,8 +874,9 @@ def test_masked_constraint_model( masked_constraint_model.solve( solver, io_api=io_api, explicit_coordinate_names=explicit_coordinate_names ) - assert (np.isclose(masked_constraint_model.solution.y[:-2], 10, rtol=1e-4)).all() - assert (np.isclose(masked_constraint_model.solution.y[-2:], 5, rtol=2e-4)).all() + tol = GPU_SOL_TOL if solver in gpu_solvers else CPU_SOL_TOL + assert (np.isclose(masked_constraint_model.solution.y[:-2], 10, rtol=tol)).all() + assert (np.isclose(masked_constraint_model.solution.y[-2:], 5, rtol=tol)).all() @pytest.mark.parametrize("solver,io_api,explicit_coordinate_names", params) @@ -954,10 +964,8 @@ def test_model_resolve( assert status == "ok" # x = -0.75, y = 3.0 - if solver == "cupdlpx": # this solver has low resolution - assert np.isclose(model.objective.value or 0, 5.25, rtol=1e-3) - else: - assert np.isclose(model.objective.value or 0, 5.25) + tol = GPU_SOL_TOL if solver in gpu_solvers else CPU_SOL_TOL + assert np.isclose(model.objective.value or 0, 5.25, rtol=tol) @pytest.mark.parametrize(