diff --git a/linopy/constraints.py b/linopy/constraints.py index d3ebef19..9b91fe25 100644 --- a/linopy/constraints.py +++ b/linopy/constraints.py @@ -603,7 +603,7 @@ def mask_func(data: pd.DataFrame) -> pd.Series: check_has_nulls(df, name=f"{self.type} {self.name}") return df - def to_polars(self) -> pl.DataFrame: + def to_polars(self, joined=True) -> pl.DataFrame: """ Convert the constraint to a polars DataFrame. @@ -631,26 +631,35 @@ def to_polars(self) -> pl.DataFrame: mask = labels_flat != -1 labels_masked = labels_flat[mask] rhs_flat = np.broadcast_to(ds["rhs"].values, ds["labels"].shape).reshape(-1) - sign_values = ds["sign"].values - sign_flat = np.broadcast_to(sign_values, ds["labels"].shape).reshape(-1) - all_same_sign = len(sign_flat) > 0 and ( - sign_flat[0] == sign_flat[-1] and (sign_flat[0] == sign_flat).all() + + short = pl.DataFrame( + { + "labels": labels_masked, + "rhs": rhs_flat[mask], + } ) - short_data: dict = { - "labels": labels_masked, - "rhs": rhs_flat[mask], - } - if all_same_sign: - short = pl.DataFrame(short_data).with_columns( - pl.lit(sign_flat[0]).cast(pl.Enum(["=", "<=", ">="])).alias("sign") - ) - else: - short_data["sign"] = pl.Series( - "sign", sign_flat[mask], dtype=pl.Enum(["=", "<=", ">="]) + all_same_sign = ( + len(sign_values) > 0 + and ((x := np.ravel(sign_values))[0] == x[-1]) + and (x[0] == sign_values).all() + ) + + short = short.with_columns( + pl.lit(np.ravel(sign_values)[0]) + .cast(pl.Enum(["=", "<=", ">="])) + .alias("sign") + if all_same_sign + else pl.Series( + "sign", + np.ravel(np.broadcast_to(sign_values, ds["labels"].shape)), + dtype=pl.Enum(["=", "<=", ">="]), ) - short = pl.DataFrame(short_data) + ) + + if not joined: + return long[["labels", "coeffs", "vars"]], short[["labels", "sign", "rhs"]] df = long.join(short, on="labels", how="inner") return df[["labels", "coeffs", "vars", "sign", "rhs"]] diff --git a/linopy/io.py b/linopy/io.py index 54090e87..a33148bb 100644 --- a/linopy/io.py +++ b/linopy/io.py @@ -902,6 +902,162 @@ def to_cupdlpx(m: Model, explicit_coordinate_names: bool = False) -> cupdlpxMode return cu_model +def to_poi( + m: Model, poi_model: Any, slice_size: int = 2_000_000 +) -> tuple[np.ndarray, np.ndarray]: + """ + Transfer a linopy model into an existing pyoptinterface model. + + Variables and constraints are added one at a time following the pyoptinterface + API. Returns two numpy arrays for result retrieval: + - ``vars_to_poi``: linopy variable label -> POI variable index + - ``cons_to_poi``: linopy constraint label -> POI constraint index + + Notes + ----- + Heavily inspired by pyoframe: https://github.com/Bravos-Power/pyoframe/blob/cf0cc39bd54e6fea0c95b0c36cd43fc0633527e1/src/pyoframe/_core.py#L1778-L1925 + + Parameters + ---------- + m : linopy.Model + poi_model : pyoptinterface model instance + An already-constructed POI model (e.g. ``pyoptinterface.highs.Model()``). + """ + import pyoptinterface as poi + + num_vars = m._xCounter or int( + max(v.data["labels"].max() for _, v in m.variables.items()) + ) + num_constrs = m._cCounter or int( + max(c.data["labels"].max() for _, c in m.constraints.items()) + ) + + # --- Variables --- + # Build a direct lookup array: linopy label -> POI variable index (O(1) per lookup) + vars_to_poi = np.full(num_vars + 1, -1, dtype=np.int32) + + for name, var in m.variables.items(): + df = var.to_polars().with_columns( + pl.col("lower").fill_nan(-np.inf).cast(pl.Float64), + pl.col("upper").fill_nan(np.inf).cast(pl.Float64), + ) # columns: lower, upper, labels + if name in m.binaries: + domain = poi.VariableDomain.Binary + elif name in m.integers: + domain = poi.VariableDomain.Integer + else: + domain = poi.VariableDomain.Continuous + + vars_to_poi[df["labels"].to_numpy()] = [ + poi_model.add_variable(domain, lower, upper, f"V{label}").index + for lower, upper, label in df.select("lower", "upper", "labels").iter_rows() + ] + + # --- Constraints --- + sense_map = { + "<=": poi.ConstraintSense.LessEqual, + ">=": poi.ConstraintSense.GreaterEqual, + "=": poi.ConstraintSense.Equal, + } + + cons_to_poi = np.full(num_constrs + 1, -1, dtype=np.int32) + + for _, con in m.constraints.items(): + for con_slice in con.iterate_slices(slice_size): + long, short = con_slice.to_polars( + joined=False + ) # columns: labels, coeffs, vars, sign, rhs + if long.is_empty(): + continue + + poi_vars = vars_to_poi[long["vars"].to_numpy()].tolist() + coeffs = long["coeffs"].to_list() + + # Compute group boundaries (first row index of each label group) + df_unique = ( + long.lazy() + .with_row_index() + .filter(pl.col("labels").is_first_distinct()) + .select("index", "labels") + .join(short.lazy(), on="labels", how="inner") + .with_columns(name=pl.lit("C") + pl.col("labels").cast(pl.String)) + .collect() + ) + + split = df_unique["index"].to_list() + [long.height] + rhs_list = df_unique["rhs"].to_list() + sign_list = df_unique["sign"].to_list() + labels = df_unique["labels"].to_numpy() + names = df_unique["name"].to_list() + + cons_to_poi[labels] = [ + poi_model.add_linear_constraint( + poi.ScalarAffineFunction(coeffs[s0:s1], poi_vars[s0:s1]), + sense_map[sign], + rhs, + name, # according to pyoframe it'd be to pass name="C" for all for gurobi + ).index + for s0, s1, rhs, sign, name in zip( + split[:-1], split[1:], rhs_list, sign_list, names + ) + ] + + # --- Objective --- + obj_df = m.objective.to_polars() + const = float(obj_df["const"].sum()) if "const" in obj_df.columns else 0.0 + obj_sense = ( + poi.ObjectiveSense.Minimize + if m.objective.sense == "min" + else poi.ObjectiveSense.Maximize + ) + if m.is_quadratic: + quad = obj_df.filter(pl.col("vars2") != -1) + obj_expr = poi.ScalarQuadraticFunction( + quad["coeffs"].to_numpy(), + vars_to_poi[quad["vars1"].to_numpy()], + vars_to_poi[quad["vars2"].to_numpy()], + ) + lin = obj_df.filter(pl.col("vars2") == -1) + if len(lin): + obj_expr = obj_expr + poi.ScalarAffineFunction.from_numpy( + lin["coeffs"].to_numpy(), vars_to_poi[lin["vars1"].to_numpy()], const + ) + poi_model.set_objective(obj_expr, obj_sense) + else: + obj_expr = poi.ScalarAffineFunction.from_numpy( + obj_df["coeffs"].to_numpy(), vars_to_poi[obj_df["vars"].to_numpy()], const + ) + poi_model.set_objective(obj_expr, obj_sense) + + # --- SOS constraints --- + sos_type_map = {1: poi.SOSType.SOS1, 2: poi.SOSType.SOS2} + for var_name in m.variables.sos: + var = m.variables.sos[var_name] + sos_type = sos_type_map[var.attrs[SOS_TYPE_ATTR]] # type: ignore[index] + sos_dim: str = var.attrs[SOS_DIM_ATTR] # type: ignore[assignment] + other_dims = [dim for dim in var.labels.dims if dim != sos_dim] + + def add_sos( + labels: xr.DataArray, sos_type: Any = sos_type, sos_dim: str = sos_dim + ) -> None: + labels = labels.squeeze() + sos_labels = labels.values.flatten() + sos_weights = labels.coords[sos_dim].values.astype(float).tolist() + poi_vars_sos = [ + poi.VariableIndex(idx) for idx in vars_to_poi[sos_labels].tolist() + ] + poi_model.add_sos_constraint(poi_vars_sos, sos_type, sos_weights) + + if not other_dims: + add_sos(var.labels) + else: + stacked = var.labels.stack(_sos_group=other_dims) + for _, s in stacked.groupby("_sos_group"): + add_sos(s.unstack("_sos_group")) + + return vars_to_poi, cons_to_poi + + 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 049093de..9215411a 100644 --- a/linopy/model.py +++ b/linopy/model.py @@ -60,6 +60,7 @@ to_highspy, to_mosek, to_netcdf, + to_poi, ) from linopy.matrices import MatrixAccessor from linopy.objective import Objective @@ -77,6 +78,7 @@ from linopy.solvers import ( IO_APIS, available_solvers, + poi_available, ) from linopy.sos_reformulation import ( reformulate_sos_constraints, @@ -1333,6 +1335,11 @@ def solve( f"Keyword argument `io_api` has to be one of {IO_APIS} or None" ) + if io_api == "poi" and not poi_available: + raise ImportError( + "pyoptinterface is not installed. Install it to use io_api='poi'." + ) + if remote is not None: if isinstance(remote, OetcHandler): solved = remote.solve_on_oetc(self) @@ -1382,17 +1389,18 @@ def solve( ) logger.info(f"Solver options:\n{options_string}") - if problem_fn is None: - problem_fn = self.get_problem_file(io_api=io_api) - if solution_fn is None: - if ( - solver_supports(solver_name, SolverFeature.SOLUTION_FILE_NOT_NEEDED) - and not keep_files - ): - # these (solver, keep_files=False) combos do not need a solution file - solution_fn = None - else: - solution_fn = self.get_solution_file() + if io_api != "poi": + if problem_fn is None: + problem_fn = self.get_problem_file(io_api=io_api) + if solution_fn is None: + if ( + solver_supports(solver_name, SolverFeature.SOLUTION_FILE_NOT_NEEDED) + and not keep_files + ): + # these (solver, keep_files=False) combos do not need a solution file + solution_fn = None + else: + solution_fn = self.get_solution_file() if sanitize_zeros: self.constraints.sanitize_zeros() @@ -1407,6 +1415,9 @@ def solve( f"Solver {solver_name} does not support quadratic problems." ) + if io_api == "poi" and not solver_supports(solver_name, SolverFeature.POI_API): + raise ValueError(f"Solver {solver_name} does not support io_api='poi'.") + if reformulate_sos not in (True, False, "auto"): raise ValueError( f"Invalid value for reformulate_sos: {reformulate_sos!r}. " @@ -1436,7 +1447,14 @@ def solve( solver = solver_class( **solver_options, ) - if io_api == "direct": + if io_api == "poi": + # no problem file written; model transferred via pyoptinterface + result = solver.solve_problem_from_poi( + model=self, + log_fn=to_path(log_fn), + env=env, + ) + elif io_api == "direct": # no problem file written and direct model is set for solver result = solver.solve_problem_from_model( model=self, @@ -1804,3 +1822,5 @@ def reset_solution(self) -> None: to_cupdlpx = to_cupdlpx to_block_files = to_block_files + + to_poi = to_poi diff --git a/linopy/poi_model.py b/linopy/poi_model.py new file mode 100644 index 00000000..5e3d3807 --- /dev/null +++ b/linopy/poi_model.py @@ -0,0 +1,216 @@ +""" +PoiModel: A wrapper that holds a pyoptinterface model alongside a linopy model, +with mappings between the two, and supports incremental updates via model diffing. +""" + +from __future__ import annotations + +from dataclasses import dataclass, field +from typing import TYPE_CHECKING, Any + +import numpy as np +import polars as pl + +if TYPE_CHECKING: + from linopy.model import Model + + +@dataclass +class PoiModel: + """ + Wraps a pyoptinterface model built from a linopy model. + + Attributes + ---------- + poi : Any + The pyoptinterface model instance. + linopy_model : Model + The linopy model used to build the POI model. + vars_to_poi : np.ndarray + Array of shape (num_vars+1,) mapping linopy variable labels to POI variable indices. + cons_to_poi : np.ndarray + Array of shape (num_constrs+1,) mapping linopy constraint labels to POI constraint indices. + """ + + poi: Any + linopy_model: Model + vars_to_poi: np.ndarray = field(repr=False) + cons_to_poi: np.ndarray = field(repr=False) + + @classmethod + def from_linopy( + cls, m: Model, poi_model: Any, slice_size: int = 2_000_000 + ) -> PoiModel: + """ + Build a PoiModel from a linopy model and an existing pyoptinterface model instance. + + Parameters + ---------- + m : linopy.Model + The linopy model to transfer. + poi_model : Any + An already-constructed POI model (e.g. ``pyoptinterface.highs.Model()``). + slice_size : int, optional + Number of constraint rows to process at once. + + Returns + ------- + PoiModel + """ + from linopy.io import to_poi + + vars_to_poi, cons_to_poi = to_poi(m, poi_model, slice_size=slice_size) + return cls( + poi=poi_model, + linopy_model=m, + vars_to_poi=vars_to_poi, + cons_to_poi=cons_to_poi, + ) + + def update_linopy(self, new_m: Model) -> None: + """ + Update the POI model in-place to reflect a new linopy model. + + Uses pyoptinterface's ``set_normalized_rhs`` and ``set_normalized_coeff`` + to apply only the differences between the stored linopy model and ``new_m``. + + Assumptions (asserted): + - Same set of constraint names. + - For each constraint, labels, variable references, and signs are identical. + - Variable labels and names are identical. + + Parameters + ---------- + new_m : linopy.Model + The updated linopy model. Must have the same structure as the original. + """ + import pyoptinterface as poi + + def make_con(idx: int) -> poi.ConstraintIndex: + return poi.ConstraintIndex(poi.ConstraintType.Linear, idx) + + old_m = self.linopy_model + + # --- Validate variable structure --- + old_vars = set(old_m.variables) + new_vars = set(new_m.variables) + assert old_vars == new_vars, ( + f"Variable sets differ. Added: {new_vars - old_vars}, " + f"Removed: {old_vars - new_vars}" + ) + + # --- Validate constraint structure --- + old_cons = set(old_m.constraints) + new_cons = set(new_m.constraints) + assert old_cons == new_cons, ( + f"Constraint sets differ. Added: {new_cons - old_cons}, " + f"Removed: {old_cons - new_cons}" + ) + + vars_to_poi = self.vars_to_poi + cons_to_poi = self.cons_to_poi + + # --- Update variable bounds --- + for vn in old_m.variables: + d1 = old_m.variables[vn].data + d2 = new_m.variables[vn].data + + assert (d1["labels"].values == d2["labels"].values).all(), ( + f"Variable '{vn}': labels differ" + ) + + lower_diff = np.ravel(d1["lower"].values != d2["lower"].values) + upper_diff = np.ravel(d1["upper"].values != d2["upper"].values) + both_diff = lower_diff & upper_diff + only_lower = lower_diff & ~upper_diff + only_upper = upper_diff & ~lower_diff + + labels = np.ravel(d2["labels"].values) + lower = np.ravel(d2["lower"].fillna(-np.inf).values) + upper = np.ravel(d2["upper"].fillna(np.inf).values) + + print(both_diff.sum(), only_lower.sum(), only_upper.sum()) + + for var_idx, lb, ub in zip( + vars_to_poi[labels[both_diff]], + lower[both_diff], + upper[both_diff], + ): + self.poi.set_variable_bounds(poi.VariableIndex(var_idx), lb, ub) + for var_idx, lb in zip( + vars_to_poi[labels[only_lower]], + lower[only_lower], + ): + self.poi.set_variable_lower_bound(poi.VariableIndex(var_idx), lb) + for var_idx, ub in zip( + vars_to_poi[labels[only_upper]], + upper[only_upper], + ): + self.poi.set_variable_upper_bound(poi.VariableIndex(var_idx), ub) + + # --- Update constraints --- + for cn in old_m.constraints: + c1 = old_m.constraints[cn].data + c2 = new_m.constraints[cn].data + + c1_labels = np.ravel(c1["labels"]) + c2_labels = np.ravel(c2["labels"]) + c1_vars = np.ravel(c1["vars"]) + c2_vars = np.ravel(c2["vars"]) + c1_sign = np.ravel(c1["sign"]) + c2_sign = np.ravel(c2["sign"]) + + assert (c1_labels == c2_labels).all(), f"Constraint '{cn}': labels differ" + assert (c1_vars == c2_vars).all(), ( + f"Constraint '{cn}': variable references differ" + ) + assert (c1_sign == c2_sign).all(), f"Constraint '{cn}': signs differ" + + # Update RHS where changed + c1_rhs = np.ravel(c1["rhs"]) + c2_rhs = np.ravel(c2["rhs"]) + rhs_diff = c1_rhs != c2_rhs + if rhs_diff.any(): + poi_cons = cons_to_poi[c2_labels[rhs_diff]] + new_rhs = c2_rhs[rhs_diff] + for con_idx, rhs_val in zip(poi_cons, new_rhs): + self.poi.set_normalized_rhs(make_con(con_idx), rhs_val) + + # Update coefficients where changed + c1_coeffs = np.ravel(c1["coeffs"]) + c2_coeffs = np.ravel(c2["coeffs"]) + # Treat NaN == NaN as equal (no change) + both_nan = np.isnan(c1_coeffs) & np.isnan(c2_coeffs) + coeffs_diff = (c1_coeffs != c2_coeffs) & ~both_nan + if coeffs_diff.any(): + broadcast_labels = np.ravel( + np.broadcast_to( + c1["labels"].values[..., np.newaxis], c1["coeffs"].shape + ) + ) + changed_labels = broadcast_labels[coeffs_diff] + changed_vars = c1_vars[coeffs_diff] + changed_coeffs = c2_coeffs[coeffs_diff] + + # Group by (con, var) and sum coefficients + updates = ( + pl.DataFrame( + { + "con": cons_to_poi[changed_labels], + "var": vars_to_poi[changed_vars], + "coeff": changed_coeffs, + } + ) + .group_by("con", "var") + .agg(pl.col("coeff").sum()) + ) + + for con_idx, var_idx, coeff_val in updates.iter_rows(): + self.poi.set_normalized_coefficient( + make_con(con_idx), + poi.VariableIndex(var_idx), + coeff_val, + ) + + # Update the stored linopy model reference + self.linopy_model = new_m diff --git a/linopy/solver_capabilities.py b/linopy/solver_capabilities.py index f0507317..0d7085ef 100644 --- a/linopy/solver_capabilities.py +++ b/linopy/solver_capabilities.py @@ -52,6 +52,7 @@ class SolverFeature(Enum): # Solver-specific SOLVER_ATTRIBUTE_ACCESS = auto() # Direct access to solver variable attributes + POI_API = auto() # Solve via pyoptinterface (poi) without writing files @dataclass(frozen=True) @@ -81,6 +82,7 @@ def supports(self, feature: SolverFeature) -> bool: SolverFeature.INTEGER_VARIABLES, SolverFeature.QUADRATIC_OBJECTIVE, SolverFeature.DIRECT_API, + SolverFeature.POI_API, SolverFeature.LP_FILE_NAMES, SolverFeature.READ_MODEL_FROM_FILE, SolverFeature.SOLUTION_FILE_NOT_NEEDED, @@ -98,6 +100,7 @@ def supports(self, feature: SolverFeature) -> bool: SolverFeature.INTEGER_VARIABLES, SolverFeature.QUADRATIC_OBJECTIVE, SolverFeature.DIRECT_API, + SolverFeature.POI_API, SolverFeature.LP_FILE_NAMES, SolverFeature.READ_MODEL_FROM_FILE, SolverFeature.SOLUTION_FILE_NOT_NEEDED, @@ -204,6 +207,7 @@ def supports(self, feature: SolverFeature) -> bool: SolverFeature.INTEGER_VARIABLES, SolverFeature.QUADRATIC_OBJECTIVE, SolverFeature.DIRECT_API, + SolverFeature.POI_API, SolverFeature.LP_FILE_NAMES, SolverFeature.READ_MODEL_FROM_FILE, SolverFeature.SOLUTION_FILE_NOT_NEEDED, @@ -217,6 +221,7 @@ def supports(self, feature: SolverFeature) -> bool: { SolverFeature.INTEGER_VARIABLES, SolverFeature.QUADRATIC_OBJECTIVE, + SolverFeature.POI_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..6cd856bf 100644 --- a/linopy/solvers.py +++ b/linopy/solvers.py @@ -52,9 +52,10 @@ ) FILE_IO_APIS = ["lp", "lp-polars", "mps"] -IO_APIS = FILE_IO_APIS + ["direct"] +IO_APIS = FILE_IO_APIS + ["direct", "poi"] available_solvers = [] +poi_available = False which = "where" if os.name == "nt" else "which" @@ -219,6 +220,10 @@ class xpress_Namespaces: # type: ignore[no-redef] except ImportError: pass +with contextlib.suppress(ModuleNotFoundError): + import pyoptinterface as poi + + poi_available = True quadratic_solvers = [s for s in QUADRATIC_SOLVERS if s in available_solvers] logger = logging.getLogger(__name__) @@ -297,6 +302,131 @@ def maybe_adjust_objective_sign( return solution +class POIMixin: + """ + Mixin for solvers that support solving via pyoptinterface. + + Subclasses must implement `_make_poi_model()` to return a configured + POI model instance (after calling autoload_library if needed). + """ + + termination_map = { + "OPTIMIZE_NOT_CALLED": "unknown", + "OPTIMAL": "optimal", + "INFEASIBLE": "infeasible", + "DUAL_INFEASIBLE": "unbounded", + "LOCALLY_SOLVED": "suboptimal", + "LOCALLY_INFEASIBLE": "infeasible", + "INFEASIBLE_OR_UNBOUNDED": "infeasible_or_unbounded", + "ALMOST_OPTIMAL": "suboptimal", + "SLOW_PROGRESS": "other", + "NUMERICAL_ERROR": "internal_solver_error", + "SOLUTION_LIMIT": "terminated_by_limit", + "ITERATION_LIMIT": "iteration_limit", + "TIME_LIMIT": "time_limit", + "NODE_LIMIT": "terminated_by_limit", + "OTHER_LIMIT": "terminated_by_limit", + "INTERRUPTED": "user_interrupt", + "OTHER_ERROR": "internal_solver_error", + } + + @abstractmethod + def _make_poi_model(self, log_fn: Path | None = None, env: Any = None) -> Any: + """Create and return a solver-specific POI model instance.""" + ... + + @abstractmethod + def safe_get_solution( + self, status: Status, func: Callable[[], Solution] + ) -> Solution: + """ + Get solution from function call, handling unknown status. + + Provided by Solver base class; declared here for type checking. + """ + ... + + def solve_problem_from_poi( + self, + model: Model, + log_fn: Path | None = None, + env: Any = None, + ) -> Result: + """ + Solve a linopy model via pyoptinterface without writing any files. + + Parameters + ---------- + model : linopy.Model + log_fn : Path, optional + Path to the log file. + env : optional + Solver-specific environment object (e.g. pyoptinterface gurobi.Env). + + Returns + ------- + Result + """ + if not poi_available: + raise ModuleNotFoundError( + "pyoptinterface is not installed. Please install it to use POI-based solvers." + ) + + poi_model = self._make_poi_model(log_fn=log_fn, env=env) + + vars_to_poi, cons_to_poi = linopy.io.to_poi(model, poi_model) + + poi_model.optimize() + + termination_status = poi_model.get_model_attribute( + poi.ModelAttribute.TerminationStatus + ) + termination_condition = POIMixin.termination_map.get( + termination_status.name, "other" + ) + status = Status.from_termination_condition(termination_condition) + status.legacy_status = termination_status + + def get_solver_solution() -> Solution: + objective = poi_model.get_model_attribute(poi.ModelAttribute.ObjectiveValue) + + # vars_to_poi / cons_to_poi: linopy label -> POI index (int32), -1 = missing + var_labels = np.nonzero(vars_to_poi != -1)[0] + var_values = np.array( + [ + poi_model.get_variable_attribute( + poi.VariableIndex(poi_var), + poi.VariableAttribute.Value, + ) + for poi_var in vars_to_poi[var_labels] + ], + dtype=float, + ) + primal = pd.Series(var_values, index=var_labels) + + try: + con_labels = np.nonzero(cons_to_poi != -1)[0] + con_values = np.array( + [ + poi_model.get_constraint_attribute( + poi.ConstraintIndex(poi.ConstraintType.Linear, poi_con), + poi.ConstraintAttribute.Dual, + ) + for poi_con in cons_to_poi[con_labels] + ], + dtype=float, + ) + dual = pd.Series(con_values, index=con_labels) + except Exception: + logger.warning("Dual values couldn't be parsed") + dual = pd.Series(dtype=float) + + return Solution(primal, dual, objective) + + solution = self.safe_get_solution(status=status, func=get_solver_solution) + return Result(status, solution, poi_model) + + class Solver(ABC, Generic[EnvType]): """ Abstract base class for solving a given linear problem. @@ -780,7 +910,7 @@ def get_solver_solution() -> Solution: return Result(status, solution) -class Highs(Solver[None]): +class Highs(Solver[None], POIMixin): """ Solver subclass for the HiGHS solver. HiGHS must be installed for usage. Find the documentation at https://highs.dev/. @@ -806,6 +936,23 @@ def __init__( ) -> None: super().__init__(**solver_options) + def _make_poi_model(self, log_fn: Path | None = None, env: Any = None) -> Any: + import pyoptinterface.highs as highs + + m = highs.Model() + for k, v in self.solver_options.items(): + if isinstance(v, bool): + m.set_raw_option_bool(k, v) + elif isinstance(v, int): + m.set_raw_option_int(k, v) + elif isinstance(v, float): + m.set_raw_option_double(k, v) + else: + m.set_raw_option_string(k, v) + if log_fn is not None: + m.set_raw_option_string("log_file", path_to_string(log_fn)) + return m + def solve_problem_from_model( self, model: Model, @@ -1031,7 +1178,7 @@ def get_solver_solution() -> Solution: return Result(status, solution, h) -class Gurobi(Solver["gurobipy.Env | dict[str, Any] | None"]): +class Gurobi(Solver["gurobipy.Env | dict[str, Any] | None"], POIMixin): """ Solver subclass for the gurobi solver. @@ -1047,6 +1194,17 @@ def __init__( ) -> None: super().__init__(**solver_options) + def _make_poi_model(self, log_fn: Path | None = None, env: Any = None) -> Any: + import pyoptinterface.gurobi as gurobi + + gurobi.autoload_library() + m = gurobi.Model(env) if env is not None else gurobi.Model() + for k, v in self.solver_options.items(): + m.set_raw_parameter(k, v) + if log_fn is not None: + m.set_raw_parameter("LogFile", path_to_string(log_fn)) + return m + def solve_problem_from_model( self, model: Model, @@ -1942,7 +2100,7 @@ def get_solver_solution() -> Solution: mosek_bas_re = re.compile(r" (XL|XU)\s+([^ \t]+)\s+([^ \t]+)| (LL|UL|BS)\s+([^ \t]+)") -class Mosek(Solver[None]): +class Mosek(Solver[None], POIMixin): """ Solver subclass for the Mosek solver. @@ -1968,6 +2126,16 @@ def __init__( ) -> None: super().__init__(**solver_options) + def _make_poi_model(self, log_fn: Path | None = None, env: Any = None) -> Any: + import pyoptinterface.mosek as mosek + + m = mosek.Model() + for k, v in self.solver_options.items(): + m.set_raw_parameter(k, v) + if log_fn is not None: + m.set_raw_parameter("MSK_SPAR_LOG_FILE_PATH", path_to_string(log_fn)) + return m + def solve_problem_from_model( self, model: Model, @@ -2281,7 +2449,7 @@ def get_solver_solution() -> Solution: return Result(status, solution) -class COPT(Solver[None]): +class COPT(Solver[None], POIMixin): """ Solver subclass for the COPT solver. @@ -2302,6 +2470,16 @@ def __init( ) -> None: super().__init__(**solver_options) + def _make_poi_model(self, log_fn: Path | None = None, env: Any = None) -> Any: + import pyoptinterface.copt as copt + + m = copt.Model() + for k, v in self.solver_options.items(): + m.set_raw_parameter(k, v) + # COPT log file requires Model.setLogFile() on the native object, + # which is not accessible via the POI interface — log_fn is ignored here. + return m + def solve_problem_from_model( self, model: Model, diff --git a/pyproject.toml b/pyproject.toml index aaac2cf1..3deb63b7 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -89,6 +89,10 @@ solvers = [ "knitro>=15.1.0", # "cupdlpx>=0.1.2", pip package currently unstable ] +poi = [ + "pyoptinterface>=0.5.1", + "highsbox>=1.13.1", +] [tool.setuptools.packages.find] include = ["linopy"]