feat: add pyoptinterface io_api

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"]]

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.

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