diff --git a/linopy/constraints.py b/linopy/constraints.py index d3ebef19..3190d19b 100644 --- a/linopy/constraints.py +++ b/linopy/constraints.py @@ -1091,6 +1091,34 @@ def flat(self) -> pd.DataFrame: df["key"] = df.labels.map(map_labels) return df + def to_polars(self) -> pl.DataFrame: + """ + Convert all constraints to a single polars DataFrame. + + The resulting dataframe is a long format with columns + `labels`, `coeffs`, `vars`, `rhs`, `sign`, `key`. + """ + dfs = [self[k].to_polars() for k in self] + if not dfs: + return pl.DataFrame( + { + "labels": pl.Series([], dtype=pl.Int64), + "coeffs": pl.Series([], dtype=pl.Float64), + "vars": pl.Series([], dtype=pl.Int64), + "sign": pl.Series([], dtype=pl.String), + "rhs": pl.Series([], dtype=pl.Float64), + "key": pl.Series([], dtype=pl.Int64), + } + ) + + df = pl.concat(dfs, how="vertical_relaxed") + labels = ( + df.select("labels") + .unique(maintain_order=True) + .with_row_index(name="key", offset=0) + ) + return df.join(labels, on="labels", how="left") + def to_matrix(self, filter_missings: bool = True) -> scipy.sparse.csc_matrix: """ Construct a constraint matrix in sparse format. @@ -1098,25 +1126,88 @@ def to_matrix(self, filter_missings: bool = True) -> scipy.sparse.csc_matrix: Missing values, i.e. -1 in labels and vars, are ignored filtered out. """ - # TODO: rename "filter_missings" to "~labels_as_coordinates" - cons = self.flat - if not len(self): raise ValueError("No constraints available to convert to matrix.") - if filter_missings: - vars = self.model.variables.flat - shape = (cons.key.max() + 1, vars.key.max() + 1) - cons["vars"] = cons.vars.map(vars.set_index("labels").key) - return scipy.sparse.csc_matrix( - (cons.coeffs, (cons.key, cons.vars)), shape=shape - ) + def _build_dense_key_map( + arrays: list[np.ndarray], total_size: int + ) -> tuple[np.ndarray, int]: + mapping = np.full(total_size, -1, dtype=np.int64) + next_key = 0 + for labels in arrays: + labels = labels[labels != -1] + if labels.size: + # Keep first-seen order while de-duplicating labels. + first_idx = np.unique(labels, return_index=True)[1] + labels = labels[np.sort(first_idx)] + n = labels.size + if n: + mapping[labels] = np.arange(next_key, next_key + n) + next_key += n + return mapping, next_key + + # Build sparse triplets directly from NumPy arrays to avoid dataframe overhead. + row_parts: list[np.ndarray] = [] + col_parts: list[np.ndarray] = [] + data_parts: list[np.ndarray] = [] + + for _, constraint in self.items(): + labels = constraint.labels.values.reshape(-1) + vars_arr = constraint.vars.values + coeffs_arr = constraint.coeffs.values + + term_axis = constraint.vars.get_axis_num(constraint.term_dim) + if term_axis != vars_arr.ndim - 1: + vars_arr = np.moveaxis(vars_arr, term_axis, -1) + coeffs_arr = np.moveaxis(coeffs_arr, term_axis, -1) + + nterm = vars_arr.shape[-1] + row = np.repeat(labels, nterm) + col = vars_arr.reshape(-1) + data = coeffs_arr.reshape(-1) + + mask = (row != -1) & (col != -1) & (data != 0) + if mask.any(): + row_parts.append(row[mask]) + col_parts.append(col[mask]) + data_parts.append(data[mask]) + + if row_parts: + row = np.concatenate(row_parts) + col = np.concatenate(col_parts) + data = np.concatenate(data_parts) else: - shape = self.model.shape - return scipy.sparse.csc_matrix( - (cons.coeffs, (cons.labels, cons.vars)), shape=shape + row = np.array([], dtype=np.int64) + col = np.array([], dtype=np.int64) + data = np.array([], dtype=float) + + if filter_missings: + # Keep the filtered matrix row space aligned with ``self.flat``: + # only constraints that still have at least one active coefficient + # get a key when ``filter_missings=True``. + cons_map, next_con_key = _build_dense_key_map([row], self.model._cCounter) + vars_map, next_var_key = _build_dense_key_map( + [ + var.labels.values.reshape(-1) + for _, var in self.model.variables.items() + ], + self.model._xCounter, ) + shape = (next_con_key, next_var_key) + + row = cons_map[row] + col = vars_map[col] + keep = (row != -1) & (col != -1) + if not keep.all(): + row = row[keep] + col = col[keep] + data = data[keep] + return scipy.sparse.csc_matrix((data, (row, col)), shape=shape) + + shape = self.model.shape + return scipy.sparse.csc_matrix((data, (row, col)), shape=shape) + def reset_dual(self) -> None: """ Reset the stored solution of variables. diff --git a/linopy/io.py b/linopy/io.py index 54090e87..e1314cb4 100644 --- a/linopy/io.py +++ b/linopy/io.py @@ -5,6 +5,7 @@ from __future__ import annotations +import contextlib import logging import shutil import time @@ -836,6 +837,271 @@ def to_highspy(m: Model, explicit_coordinate_names: bool = False) -> Highs: return h +def to_xpress( + m: Model, explicit_coordinate_names: bool = False, progress: bool | None = None +) -> Any: + """ + Export the model to xpress using native array-loading APIs. + + The model is transferred through loadLP/loadQP/loadMIQP, matching the + underlying Xpress Optimizer C API data layout. + """ + import xpress + + if progress is None: + progress = m._xCounter > 10_000 + + def _name_array(labels: np.ndarray, formatter: Callable[[Any], str]) -> np.ndarray: + flat_labels = labels.ravel() + return np.fromiter( + (formatter(label) for label in flat_labels), + dtype=object, + count=flat_labels.size, + ) + + M = m.matrices + t_start = time.perf_counter() + t_stage = t_start + problem = xpress.problem() + vlabels = M.vlabels + clabels = M.clabels + A = M.A + Q = M.Q + int32_max = np.iinfo(np.int32).max + use_int32_indices = ( + len(vlabels) <= int32_max + and len(clabels) <= int32_max + and (A is None or A.nnz <= int32_max) + and (Q is None or Q.nnz <= int32_max) + ) + index_dtype = np.int32 if use_int32_indices else np.int64 + + def _emit_progress_message(message: str) -> None: + if not progress: + return + logger.info(message) + + _emit_progress_message( + " Xpress direct IO: building model " + f"(nvars={len(vlabels)}, ncons={len(clabels)}, " + f"annz={int(A.nnz) if A is not None else 0}, " + f"qnnz={int(Q.nnz) if Q is not None else 0}, " + f"explicit_names={explicit_coordinate_names})" + ) + + def call_xpress(new_api: str, old_api: str, **kwargs: Any) -> None: + try: + getattr(problem, new_api)(**kwargs) + except AttributeError: + getattr(problem, old_api)(**kwargs) + + def add_names(namespace: int, names: list[str]) -> None: + if not names: + return + try: + problem.addNames(namespace, names, 0, len(names) - 1) + except AttributeError: + problem.addnames(namespace, names, 0, len(names) - 1) + + def _log_stage(stage: str, detail: str = "") -> None: + if not progress: + return + nonlocal t_stage + now = time.perf_counter() + msg = ( + f" Xpress direct IO: {stage} ({now - t_stage:.3f}s stage, " + f"{now - t_start:.3f}s total)" + ) + if detail: + msg = f"{msg} {detail}" + _emit_progress_message(msg) + t_stage = now + + if A is not None and A.nnz: + if A.format != "csc": + A = A.tocsc() + start = A.indptr.astype(index_dtype, copy=False) + rowind = A.indices.astype(index_dtype, copy=False) + rowcoef = A.data.astype(float, copy=False) + else: + start = None + rowind = None + rowcoef = None + + _log_stage("prepared linear constraint matrix") + + lb = np.asarray(M.lb, dtype=float) + ub = np.asarray(M.ub, dtype=float) + + lb_inf = np.isneginf(lb) + if lb_inf.any(): + lb = lb.copy() + lb[lb_inf] = -xpress.infinity + + ub_inf = np.isposinf(ub) + if ub_inf.any(): + ub = ub.copy() + ub[ub_inf] = xpress.infinity + + _log_stage("prepared variable bounds") + + if len(clabels): + sense = M.sense + rowtype = np.full(sense.shape, "E", dtype="U1") + rowtype[sense == "<"] = "L" + rowtype[sense == ">"] = "G" + rhs = np.asarray(M.b, dtype=float) + else: + rowtype = None + rhs = None + + _log_stage("prepared row senses and rhs") + + objqcol1: np.ndarray | None + objqcol2: np.ndarray | None + objqcoef: np.ndarray | None + if Q is not None and Q.nnz: + if Q.format == "coo": # codespell:ignore coo + mask = Q.row <= Q.col + objqcol1 = Q.row[mask].astype(index_dtype, copy=False) + objqcol2 = Q.col[mask].astype(index_dtype, copy=False) + objqcoef = Q.data[mask].astype(float, copy=False) + else: + Qt = triu(Q, format="coo") # codespell:ignore coo + objqcol1 = Qt.row.astype(index_dtype, copy=False) + objqcol2 = Qt.col.astype(index_dtype, copy=False) + objqcoef = Qt.data.astype(float, copy=False) + else: + objqcol1 = None + objqcol2 = None + objqcoef = None + + _log_stage("prepared quadratic objective terms") + + vtypes = M.vtypes + integer_mask = (vtypes == "B") | (vtypes == "I") + is_mip = bool(np.any(integer_mask)) + objcoef = np.asarray(M.c, dtype=float) + + if is_mip: + entind = np.flatnonzero(integer_mask).astype(index_dtype, copy=False) + coltype = vtypes[entind] + call_xpress( + "loadMIQP", + "loadmiqp", + probname="", + rowtype=rowtype, + rhs=rhs, + rng=None, + objcoef=objcoef, + start=start, + collen=None, + rowind=rowind, + rowcoef=rowcoef, + lb=lb, + ub=ub, + objqcol1=objqcol1, + objqcol2=objqcol2, + objqcoef=objqcoef, + coltype=coltype, + entind=entind, + limit=None, + settype=None, + setstart=None, + setind=None, + refval=None, + ) + elif objqcoef is not None: + call_xpress( + "loadQP", + "loadqp", + probname="", + rowtype=rowtype, + rhs=rhs, + rng=None, + objcoef=objcoef, + start=start, + collen=None, + rowind=rowind, + rowcoef=rowcoef, + lb=lb, + ub=ub, + objqcol1=objqcol1, + objqcol2=objqcol2, + objqcoef=objqcoef, + ) + else: + call_xpress( + "loadLP", + "loadlp", + probname="", + rowtype=rowtype, + rhs=rhs, + rng=None, + objcoef=objcoef, + start=start, + collen=None, + rowind=rowind, + rowcoef=rowcoef, + lb=lb, + ub=ub, + ) + + _log_stage("loaded matrix data into Xpress") + + if m.objective.sense == "max": + changed_sense = False + with contextlib.suppress(AttributeError): + problem.chgObjSense(xpress.ObjSense.MAXIMIZE) + changed_sense = True + if not changed_sense: + with contextlib.suppress(AttributeError): + problem.chgobjsense(xpress.maximize) + + _log_stage("set objective sense") + + if explicit_coordinate_names: + print_variable, print_constraint = get_printers_scalar( + m, explicit_coordinate_names=explicit_coordinate_names + ) + row_namespace = getattr(getattr(xpress, "Namespaces", None), "ROW", 1) + col_namespace = getattr(getattr(xpress, "Namespaces", None), "COLUMN", 2) + + add_names(col_namespace, _name_array(vlabels, print_variable).tolist()) + add_names(row_namespace, _name_array(clabels, print_constraint).tolist()) + + _log_stage("attached variable/constraint names") + + if m.variables.sos: + for var_name in m.variables.sos: + var = m.variables.sos[var_name] + sos_type: int = var.attrs[SOS_TYPE_ATTR] # type: ignore[assignment] + sos_dim: str = var.attrs[SOS_DIM_ATTR] # type: ignore[assignment] + + def add_sos(s: xr.DataArray, sos_type: int, sos_dim: str) -> None: + s = s.squeeze() + labels = s.values.astype(index_dtype, copy=False).flatten() + mask = labels != -1 + if not mask.any(): + return + indices = labels[mask].tolist() + weights = s.coords[sos_dim].values[mask].tolist() + problem.addSOS(indices, weights, type=sos_type) + + others = [dim for dim in var.labels.dims if dim != sos_dim] + if not others: + add_sos(var.labels, sos_type, sos_dim) + else: + stacked = var.labels.stack(_sos_group=others) + for _, s in stacked.groupby("_sos_group"): + add_sos(s.unstack("_sos_group"), sos_type, sos_dim) + _log_stage("attached SOS constraints") + + _log_stage("finished direct model build") + + return problem + + def to_cupdlpx(m: Model, explicit_coordinate_names: bool = False) -> cupdlpxModel: """ Export the model to cupdlpx. diff --git a/linopy/matrices.py b/linopy/matrices.py index a55bb0bd..2950d3d4 100644 --- a/linopy/matrices.py +++ b/linopy/matrices.py @@ -18,6 +18,7 @@ from scipy.sparse._csc import csc_matrix from linopy import expressions +from linopy.constants import FACTOR_DIM if TYPE_CHECKING: from linopy.model import Model @@ -51,11 +52,136 @@ def __init__(self, model: Model) -> None: def clean_cached_properties(self) -> None: """Clear the cache for all cached properties of an object""" - for cached_prop in ["flat_vars", "flat_cons", "sol", "dual"]: + for cached_prop in [ + "flat_vars", + "flat_cons", + "sol", + "dual", + "_variable_data", + "_constraint_data", + ]: # check existence of cached_prop without creating it if cached_prop in self.__dict__: delattr(self, cached_prop) + @cached_property + def _variable_data(self) -> tuple[ndarray, ndarray, ndarray, ndarray, ndarray]: + """Dense-by-key variable vectors and label->key map.""" + m = self._parent + + label_to_key = np.full(m._xCounter, -1, dtype=np.int64) + labels_parts: list[np.ndarray] = [] + lb_parts: list[np.ndarray] = [] + ub_parts: list[np.ndarray] = [] + vtype_parts: list[np.ndarray] = [] + next_key = 0 + + for _, variable in m.variables.items(): + labels = variable.labels.values.reshape(-1) + mask = labels != -1 + labels = labels[mask] + n = labels.size + if not n: + continue + + label_to_key[labels] = np.arange(next_key, next_key + n) + next_key += n + + lb = np.broadcast_to(variable.lower.values, variable.labels.shape).reshape( + -1 + ) + ub = np.broadcast_to(variable.upper.values, variable.labels.shape).reshape( + -1 + ) + + labels_parts.append(labels) + lb_parts.append(lb[mask]) + ub_parts.append(ub[mask]) + + if variable.attrs["binary"]: + vtype = "B" + elif variable.attrs["integer"]: + vtype = "I" + else: + vtype = "C" + vtype_parts.append(np.full(n, vtype, dtype=" tuple[ndarray, ndarray, ndarray, ndarray]: + """Dense-by-key constraint vectors and label->key map.""" + m = self._parent + + label_to_key = np.full(m._cCounter, -1, dtype=np.int64) + labels_parts: list[np.ndarray] = [] + rhs_parts: list[np.ndarray] = [] + sense_parts: list[np.ndarray] = [] + next_key = 0 + + for _, constraint in m.constraints.items(): + labels = constraint.labels.values.reshape(-1) + vars_arr = constraint.vars.values + coeffs_arr = constraint.coeffs.values + + term_axis = constraint.vars.get_axis_num(constraint.term_dim) + if term_axis != vars_arr.ndim - 1: + vars_arr = np.moveaxis(vars_arr, term_axis, -1) + coeffs_arr = np.moveaxis(coeffs_arr, term_axis, -1) + + active = ((vars_arr != -1) & (coeffs_arr != 0)).any(axis=-1).reshape(-1) + mask = (labels != -1) & active + labels = labels[mask] + n = labels.size + if not n: + continue + + label_to_key[labels] = np.arange(next_key, next_key + n) + next_key += n + + rhs = np.broadcast_to( + constraint.rhs.values, constraint.labels.shape + ).reshape(-1) + sign = np.broadcast_to( + constraint.sign.values, constraint.labels.shape + ).reshape(-1) + sign = sign[mask] + sense = np.full(sign.shape, "=", dtype="="] = ">" + + labels_parts.append(labels) + rhs_parts.append(rhs[mask]) + sense_parts.append(sense) + + if labels_parts: + return ( + np.concatenate(labels_parts), + np.concatenate(sense_parts), + np.concatenate(rhs_parts), + label_to_key, + ) + return ( + np.array([], dtype=np.int64), + np.array([], dtype=" pd.DataFrame: m = self._parent @@ -69,33 +195,20 @@ def flat_cons(self) -> pd.DataFrame: @property def vlabels(self) -> ndarray: """Vector of labels of all non-missing variables.""" - df: pd.DataFrame = self.flat_vars - return create_vector(df.key, df.labels, -1) + labels, *_ = self._variable_data + return labels @property def vtypes(self) -> ndarray: """Vector of types of all non-missing variables.""" - m = self._parent - df: pd.DataFrame = self.flat_vars - specs = [] - for name in m.variables: - if name in m.binaries: - val = "B" - elif name in m.integers: - val = "I" - else: - val = "C" - specs.append(pd.Series(val, index=m.variables[name].flat.labels)) - - ds = pd.concat(specs) - ds = df.set_index("key").labels.map(ds) - return create_vector(ds.index, ds.to_numpy(), fill_value="") + _, _, _, vtypes, _ = self._variable_data + return vtypes @property def lb(self) -> ndarray: """Vector of lower bounds of all non-missing variables.""" - df: pd.DataFrame = self.flat_vars - return create_vector(df.key, df.lower) + _, lb, _, _, _ = self._variable_data + return lb @cached_property def sol(self) -> ndarray: @@ -124,16 +237,14 @@ def dual(self) -> ndarray: @property def ub(self) -> ndarray: """Vector of upper bounds of all non-missing variables.""" - df: pd.DataFrame = self.flat_vars - return create_vector(df.key, df.upper) + _, _, ub, _, _ = self._variable_data + return ub @property def clabels(self) -> ndarray: """Vector of labels of all non-missing constraints.""" - df: pd.DataFrame = self.flat_cons - if df.empty: - return np.array([], dtype=int) - return create_vector(df.key, df.labels, fill_value=-1) + labels, _, _, _ = self._constraint_data + return labels @property def A(self) -> csc_matrix | None: @@ -141,33 +252,57 @@ def A(self) -> csc_matrix | None: m = self._parent if not len(m.constraints): return None - A: csc_matrix = m.constraints.to_matrix(filter_missings=False) - return A[self.clabels][:, self.vlabels] + return m.constraints.to_matrix(filter_missings=True) @property def sense(self) -> ndarray: """Vector of senses of all non-missing constraints.""" - df: pd.DataFrame = self.flat_cons - return create_vector(df.key, df.sign.astype(np.dtype(" ndarray: """Vector of right-hand-sides of all non-missing constraints.""" - df: pd.DataFrame = self.flat_cons - return create_vector(df.key, df.rhs) + _, _, rhs, _ = self._constraint_data + return rhs @property def c(self) -> ndarray: """Vector of objective coefficients of all non-missing variables.""" m = self._parent - ds = m.objective.flat - if isinstance(m.objective.expression, expressions.QuadraticExpression): - ds = ds[(ds.vars1 == -1) | (ds.vars2 == -1)] - ds["vars"] = ds.vars1.where(ds.vars1 != -1, ds.vars2) - - vars: pd.Series = ds.vars.map(self.flat_vars.set_index("labels").key) - shape: int = self.flat_vars.key.max() + 1 - return create_vector(vars, ds.coeffs, fill_value=0.0, shape=shape) + _, _, _, _, label_to_key = self._variable_data + nvars = len(self.vlabels) + if nvars == 0: + return np.array([], dtype=float) + + expr = m.objective.expression + + if isinstance(expr, expressions.QuadraticExpression): + vars_arr = expr.data.vars.values + coeffs = expr.data.coeffs.values.reshape(-1) + factor_axis = expr.data.vars.get_axis_num(FACTOR_DIM) + + vars1 = np.take(vars_arr, 0, axis=factor_axis).reshape(-1) + vars2 = np.take(vars_arr, 1, axis=factor_axis).reshape(-1) + mask = ((vars1 == -1) ^ (vars2 == -1)) & (coeffs != 0) + lin_vars = np.where(vars1 == -1, vars2, vars1) + labels = lin_vars[mask] + coeffs = coeffs[mask] + else: + vars_arr = expr.vars.values.reshape(-1) + coeffs = expr.coeffs.values.reshape(-1) + mask = (vars_arr != -1) & (coeffs != 0) + labels = vars_arr[mask] + coeffs = coeffs[mask] + + keys = label_to_key[labels] + valid = keys != -1 + if not np.any(valid): + return np.zeros(nvars, dtype=float) + return np.bincount(keys[valid], weights=coeffs[valid], minlength=nvars).astype( + float, + copy=False, + ) @property def Q(self) -> csc_matrix | None: diff --git a/linopy/model.py b/linopy/model.py index 049093de..60507928 100644 --- a/linopy/model.py +++ b/linopy/model.py @@ -60,6 +60,7 @@ to_highspy, to_mosek, to_netcdf, + to_xpress, ) from linopy.matrices import MatrixAccessor from linopy.objective import Objective @@ -1801,6 +1802,8 @@ def reset_solution(self) -> None: to_highspy = to_highspy + to_xpress = to_xpress + to_cupdlpx = to_cupdlpx to_block_files = to_block_files diff --git a/linopy/solver_capabilities.py b/linopy/solver_capabilities.py index f0507317..46f2cc39 100644 --- a/linopy/solver_capabilities.py +++ b/linopy/solver_capabilities.py @@ -144,6 +144,7 @@ def supports(self, feature: SolverFeature) -> bool: { SolverFeature.INTEGER_VARIABLES, SolverFeature.QUADRATIC_OBJECTIVE, + SolverFeature.DIRECT_API, SolverFeature.LP_FILE_NAMES, SolverFeature.READ_MODEL_FROM_FILE, SolverFeature.SOLUTION_FILE_NOT_NEEDED, @@ -154,6 +155,7 @@ def supports(self, feature: SolverFeature) -> bool: else { SolverFeature.INTEGER_VARIABLES, SolverFeature.QUADRATIC_OBJECTIVE, + SolverFeature.DIRECT_API, SolverFeature.LP_FILE_NAMES, SolverFeature.READ_MODEL_FROM_FILE, SolverFeature.SOLUTION_FILE_NOT_NEEDED, diff --git a/linopy/solvers.py b/linopy/solvers.py index 16c07932..cd26bf9f 100644 --- a/linopy/solvers.py +++ b/linopy/solvers.py @@ -14,6 +14,7 @@ import subprocess as sub import sys import threading +import time import warnings from abc import ABC, abstractmethod from collections import namedtuple @@ -1604,8 +1605,39 @@ def solve_problem_from_model( env: None = None, explicit_coordinate_names: bool = False, ) -> Result: - msg = "Direct API not implemented for Xpress" - raise NotImplementedError(msg) + variable_names = ( + np.asarray(model.matrices.vlabels) + if not explicit_coordinate_names + else None + ) + constraint_names = ( + np.asarray(model.matrices.clabels) + if not explicit_coordinate_names + else None + ) + + build_start = time.perf_counter() + logger.info(" Start building Xpress direct model") + m = model.to_xpress( + explicit_coordinate_names=explicit_coordinate_names, + progress=None, + ) + logger.info( + " Finished building Xpress direct model in %.3fs", + time.perf_counter() - build_start, + ) + + return self._solve( + m=m, + solution_fn=solution_fn, + log_fn=log_fn, + warmstart_fn=warmstart_fn, + basis_fn=basis_fn, + io_api="direct", + sense=model.sense, + variable_names=variable_names, + constraint_names=constraint_names, + ) def solve_problem_from_file( self, @@ -1643,6 +1675,45 @@ def solve_problem_from_file( ------- Result """ + m = xpress.problem() + + try: # Try new API first + m.readProb(path_to_string(problem_fn)) + except AttributeError: # Fallback to old API + m.read(path_to_string(problem_fn)) + + return self._solve( + m=m, + solution_fn=solution_fn, + log_fn=log_fn, + warmstart_fn=warmstart_fn, + basis_fn=basis_fn, + io_api=read_io_api_from_problem_file(problem_fn), + sense=read_sense_from_problem_file(problem_fn), + ) + + def _solve( + self, + m: Any, + solution_fn: Path | None, + log_fn: Path | None, + warmstart_fn: Path | None, + basis_fn: Path | None, + io_api: str | None, + sense: str | None, + variable_names: np.ndarray | None = None, + constraint_names: np.ndarray | None = None, + ) -> Result: + def _get_names( + namespace: int, count: int, cached: np.ndarray | None = None + ) -> np.ndarray | list[str]: + if cached is not None: + return cached + try: # Try new API first + return m.getNameList(namespace, 0, count - 1) + except AttributeError: # Fallback to old API + return m.getnamelist(namespace, 0, count - 1) + CONDITION_MAP = { xpress.SolStatus.NOTFOUND: "unknown", xpress.SolStatus.OPTIMAL: "optimal", @@ -1651,16 +1722,6 @@ def solve_problem_from_file( xpress.SolStatus.UNBOUNDED: "unbounded", } - io_api = read_io_api_from_problem_file(problem_fn) - sense = read_sense_from_problem_file(problem_fn) - - m = xpress.problem() - - try: # Try new API first - m.readProb(path_to_string(problem_fn)) - except AttributeError: # Fallback to old API - m.read(path_to_string(problem_fn)) - # Set solver options - new API uses setControl per option, old API accepts dict if self.solver_options is not None: m.setControl(self.solver_options) @@ -1677,7 +1738,12 @@ def solve_problem_from_file( except AttributeError: # Fallback to old API m.readbasis(path_to_string(warmstart_fn)) + optimize_start = time.perf_counter() + logger.info(" Start Xpress optimize()") m.optimize() + logger.info( + " Finished Xpress optimize() in %.3fs", time.perf_counter() - optimize_start + ) # if the solver is stopped (timelimit for example), postsolve the problem if m.attributes.solvestatus == xpress.enums.SolveStatus.STOPPED: @@ -1712,10 +1778,9 @@ def solve_problem_from_file( def get_solver_solution() -> Solution: objective = m.attributes.objval - try: # Try new API first - var = m.getNameList(xpress_Namespaces.COLUMN, 0, m.attributes.cols - 1) - except AttributeError: # Fallback to old API - var = m.getnamelist(xpress_Namespaces.COLUMN, 0, m.attributes.cols - 1) + var = _get_names( + xpress_Namespaces.COLUMN, m.attributes.cols, cached=variable_names + ) sol = pd.Series(m.getSolution(), index=var, dtype=float) try: @@ -1724,14 +1789,11 @@ def get_solver_solution() -> Solution: except AttributeError: # Fallback to old API _dual = m.getDual() - try: # Try new API first - constraints = m.getNameList( - xpress_Namespaces.ROW, 0, m.attributes.rows - 1 - ) - except AttributeError: # Fallback to old API - constraints = m.getnamelist( - xpress_Namespaces.ROW, 0, m.attributes.rows - 1 - ) + constraints = _get_names( + xpress_Namespaces.ROW, + m.attributes.rows, + cached=constraint_names, + ) dual = pd.Series(_dual, index=constraints, dtype=float) except (xpress.SolverError, xpress.ModelError, SystemError): logger.warning("Dual values of MILP couldn't be parsed") diff --git a/test/test_io.py b/test/test_io.py index e8ded144..eba3f52b 100644 --- a/test/test_io.py +++ b/test/test_io.py @@ -6,6 +6,8 @@ """ import pickle +import sys +import types from pathlib import Path import numpy as np @@ -15,7 +17,7 @@ import xarray as xr from linopy import LESS_EQUAL, Model, available_solvers, read_netcdf -from linopy.io import signed_number +from linopy.io import signed_number, to_xpress from linopy.testing import assert_model_equal @@ -204,6 +206,114 @@ def test_to_highspy(model: Model) -> None: model.to_highspy() +class _FakeXpressProblem: + def __init__(self) -> None: + self.calls: list[tuple[str, dict]] = [] + + def loadLP(self, **kwargs: object) -> None: + self.calls.append(("loadLP", kwargs)) + + def loadQP(self, **kwargs: object) -> None: + self.calls.append(("loadQP", kwargs)) + + def loadMIQP(self, **kwargs: object) -> None: + self.calls.append(("loadMIQP", kwargs)) + + def chgObjSense(self, objsense: object) -> None: + self.calls.append(("chgObjSense", {"objsense": objsense})) + + def addNames( + self, namespace_type: int, names: list[str], first: int, last: int + ) -> None: + self.calls.append( + ( + "addNames", + { + "type": namespace_type, + "names": names, + "first": first, + "last": last, + }, + ) + ) + + +def _install_fake_xpress(monkeypatch: pytest.MonkeyPatch) -> _FakeXpressProblem: + problem = _FakeXpressProblem() + fake_xpress = types.ModuleType("xpress") + fake_xpress.infinity = 1.0e20 + fake_xpress.ObjSense = types.SimpleNamespace(MAXIMIZE="MAX", MINIMIZE="MIN") + fake_xpress.Namespaces = types.SimpleNamespace(ROW=1, COLUMN=2) + fake_xpress.problem = lambda: problem # type: ignore[assignment] + monkeypatch.setitem(sys.modules, "xpress", fake_xpress) + return problem + + +def test_to_xpress_loadlp_and_maxsense(monkeypatch: pytest.MonkeyPatch) -> None: + fake_problem = _install_fake_xpress(monkeypatch) + + m = Model() + x = m.add_variables(lower=0, upper=10, name="x") + y = m.add_variables(lower=0, upper=5, name="y") + m.add_constraints(2 * x + y, LESS_EQUAL, 12, name="c") + m.add_objective(3 * x + y, sense="max") + + solver_model = to_xpress(m, explicit_coordinate_names=True) + + assert solver_model is fake_problem + method_names = [name for name, _ in fake_problem.calls] + assert "loadLP" in method_names + assert "chgObjSense" in method_names + assert method_names.count("addNames") == 2 + + call_map = {name: kwargs for name, kwargs in fake_problem.calls} + lp_kwargs = call_map["loadLP"] + assert lp_kwargs["start"].dtype == np.int32 + assert lp_kwargs["rowind"].dtype == np.int32 + + +def test_to_xpress_loadmiqp(monkeypatch: pytest.MonkeyPatch) -> None: + fake_problem = _install_fake_xpress(monkeypatch) + + m = Model() + x = m.add_variables(lower=0, upper=10, name="x") + y = m.add_variables(binary=True, name="y") + m.add_constraints(x + y, LESS_EQUAL, 3, name="c") + m.add_objective(x * x + 2 * y) + + to_xpress(m) + + call_map = {name: kwargs for name, kwargs in fake_problem.calls} + assert "loadMIQP" in call_map + miqp_kwargs = call_map["loadMIQP"] + assert miqp_kwargs["entind"] is not None + assert miqp_kwargs["objqcoef"] is not None + assert miqp_kwargs["entind"].dtype == np.int32 + assert miqp_kwargs["objqcol1"].dtype == np.int32 + assert miqp_kwargs["objqcol2"].dtype == np.int32 + assert "addNames" not in call_map + + +def test_to_xpress_progress_logging(monkeypatch: pytest.MonkeyPatch) -> None: + _install_fake_xpress(monkeypatch) + messages: list[str] = [] + + def _log_info(msg: str, *args: object, **kwargs: object) -> None: + messages.append(msg % args if args else msg) + + monkeypatch.setattr("linopy.io.logger.info", _log_info) + + m = Model() + x = m.add_variables(lower=0, upper=10, coords=[range(3)], name="x") + m.add_constraints(x.sum(), LESS_EQUAL, 10, name="c") + m.add_objective((2 * x).sum()) + + to_xpress(m, progress=True) + + assert any("prepared linear constraint matrix" in msg for msg in messages) + assert any("finished direct model build" in msg for msg in messages) + + def test_to_blocks(tmp_path: Path) -> None: m: Model = Model() diff --git a/test/test_matrices.py b/test/test_matrices.py index 98a73564..58ada40b 100644 --- a/test/test_matrices.py +++ b/test/test_matrices.py @@ -77,3 +77,18 @@ def test_matrices_float_c() -> None: c = m.matrices.c assert np.all(c == np.array([1.5, 1.5])) + + +def test_matrices_A_with_sparse_constraint_labels() -> None: + m = Model() + + x = m.add_variables(0, 1, coords=[range(6)], name="x") + mask = xr.DataArray([True, False, True, False, True, False], dims=["dim_0"]) + m.add_constraints(x, GREATER_EQUAL, 0, mask=mask) + m.add_objective(x.sum()) + + A_full = m.constraints.to_matrix(filter_missings=False) + expected = A_full[m.matrices.clabels][:, m.matrices.vlabels] + + assert m.matrices.A is not None + assert (m.matrices.A != expected).nnz == 0