Solvers

Solvers are thin wrappers that put the numerical backends, such as Gurobi and HiGHS, behind a common interface.

Solver

Bases: object

Abstract base class for optimization solvers.

This class defines the common interface for all solver backends in FOMTools, including MIP solvers (Gurobi, HiGHS) and SDP solvers (Clarabel). Subclasses must implement the abstract methods to interface with their specific solver.

Attributes:

Name	Type	Description
vector_var_map	dict	Mapping from LinearExpression objects to their corresponding solver variables.

Source code in fomtools/solvers/solver.py

class Solver(object):
    """
    Abstract base class for optimization solvers.

    This class defines the common interface for all solver backends in FOMTools,
    including MIP solvers (Gurobi, HiGHS) and SDP solvers (Clarabel). Subclasses
    must implement the abstract methods to interface with their specific solver.

    Attributes
    ----------
    vector_var_map : dict
        Mapping from LinearExpression objects to their corresponding solver variables.
    """

    def __init__(self, solver_params=None):
        """
        Initialize the solver with optional parameters.

        Parameters
        ----------
        solver_params : dict, optional
            Dictionary of solver-specific parameters to configure.
            Keys are parameter names and values are the settings.
        """
        self._relaxed_model = None
        self.vector_var_map = {}
        # Retain the params so the attribute always exists (and copies can carry
        # it over); subclasses such as ClarabelSolver also read it at solve time.
        self.solver_params = solver_params
        if solver_params is not None:
            self.set_solver_params(solver_params)

    # Name of the backend's verbosity option, if it has one (set by subclasses).
    VERBOSITY_KEY: str | None = None

    @classmethod
    def verbosity_params(cls, params: dict | None) -> dict | None:
        """
        The verbosity option from ``params``, if it names this backend's option.

        Used to pass only the output setting of the main solver's options on to a
        secondary solver of the same backend (e.g. the OBBT LP).
        """
        if not params or cls.VERBOSITY_KEY is None or cls.VERBOSITY_KEY not in params:
            return None
        return {cls.VERBOSITY_KEY: params[cls.VERBOSITY_KEY]}

    def set_solver_params(self, solver_params):
        """
        Apply a dictionary of solver parameters.

        Parameters
        ----------
        solver_params : dict
            Dictionary mapping parameter names to their values.
            Each key-value pair is passed to ``set_param``.
        """
        for k, v in solver_params.items():
            self.set_param(k, v)

    def set_param(self, param, value):
        """
        Set an individual solver parameter.

        Parameters
        ----------
        param : str
            Name of the solver parameter.
        value : int or float
            Value to assign to the parameter.

        Raises
        ------
        NotImplementedError
            Must be implemented by subclasses.
        """
        raise NotImplementedError

    def add_param_stack_var(self, param_stack_var, param_list):
        """
        Create a stacked variable from a list of parameter variables.

        Parameters
        ----------
        param_stack_var : LinearExpression
            The new variable to register as the stacked result.
        param_list : list[LinearExpression]
            List of existing variables to stack together.

        Raises
        ------
        NotImplementedError
            Must be implemented by subclasses.
        """
        raise NotImplementedError

    def add_constr(self, constr):
        """
        Add a single constraint to the model.

        Parameters
        ----------
        constr : Constraint
            The constraint to add to the solver model.

        Raises
        ------
        NotImplementedError
            Must be implemented by subclasses.
        """
        raise NotImplementedError

    def add_var(self, iterate, lb=-np.inf, ub=np.inf):
        """
        Add a decision variable to the model.

        Parameters
        ----------
        iterate : LinearExpression
            The variable to register with the solver.
        lb : float, default=-np.inf
            Lower bound for the variable.
        ub : float, default=np.inf
            Upper bound for the variable.

        Raises
        ------
        NotImplementedError
            Must be implemented by subclasses.
        """
        raise NotImplementedError

    def constraints_to_solver(self, constraints):
        """
        Add multiple constraints to the model.

        Parameters
        ----------
        constraints : list[Constraint]
            List of constraints to add to the solver model.

        Raises
        ------
        NotImplementedError
            Must be implemented by subclasses.
        """
        raise NotImplementedError

    def optimize(self):
        """
        Solve the optimization model.

        Raises
        ------
        NotImplementedError
            Must be implemented by subclasses.
        """
        raise NotImplementedError

    def copy(self) -> "Solver":
        """
        Create an independent copy of this solver.

        The copy can be modified without affecting the original.

        Returns
        -------
        Solver
            A new solver instance with copied state.

        Raises
        ------
        NotImplementedError
            Must be implemented by subclasses.
        """
        raise NotImplementedError

__init__

__init__(solver_params=None)

Parameters:

Name	Type	Description	Default
solver_params	dict	Dictionary of solver-specific parameters to configure. Keys are parameter names and values are the settings.	None

Source code in fomtools/solvers/solver.py

def __init__(self, solver_params=None):
    """
    Initialize the solver with optional parameters.

    Parameters
    ----------
    solver_params : dict, optional
        Dictionary of solver-specific parameters to configure.
        Keys are parameter names and values are the settings.
    """
    self._relaxed_model = None
    self.vector_var_map = {}
    # Retain the params so the attribute always exists (and copies can carry
    # it over); subclasses such as ClarabelSolver also read it at solve time.
    self.solver_params = solver_params
    if solver_params is not None:
        self.set_solver_params(solver_params)

add_constr

add_constr(constr)

Add a single constraint to the model.

Parameters:

Name	Type	Description	Default
constr	Constraint	The constraint to add to the solver model.	required

Raises:

Type	Description
NotImplementedError	Must be implemented by subclasses.

Source code in fomtools/solvers/solver.py

def add_constr(self, constr):
    """
    Add a single constraint to the model.

    Parameters
    ----------
    constr : Constraint
        The constraint to add to the solver model.

    Raises
    ------
    NotImplementedError
        Must be implemented by subclasses.
    """
    raise NotImplementedError

add_param_stack_var

add_param_stack_var(param_stack_var, param_list)

Create a stacked variable from a list of parameter variables.

Parameters:

Name	Type	Description	Default
param_stack_var	LinearExpression	The new variable to register as the stacked result.	required
param_list	list[LinearExpression]	List of existing variables to stack together.	required

Raises:

Type	Description
NotImplementedError	Must be implemented by subclasses.

Source code in fomtools/solvers/solver.py

def add_param_stack_var(self, param_stack_var, param_list):
    """
    Create a stacked variable from a list of parameter variables.

    Parameters
    ----------
    param_stack_var : LinearExpression
        The new variable to register as the stacked result.
    param_list : list[LinearExpression]
        List of existing variables to stack together.

    Raises
    ------
    NotImplementedError
        Must be implemented by subclasses.
    """
    raise NotImplementedError

add_var

add_var(iterate, lb=-inf, ub=inf)

Add a decision variable to the model.

Parameters:

Name	Type	Description	Default
iterate	LinearExpression	The variable to register with the solver.	required
lb	float	Lower bound for the variable.	-np.inf
ub	float	Upper bound for the variable.	np.inf

Raises:

Type	Description
NotImplementedError	Must be implemented by subclasses.

Source code in fomtools/solvers/solver.py

def add_var(self, iterate, lb=-np.inf, ub=np.inf):
    """
    Add a decision variable to the model.

    Parameters
    ----------
    iterate : LinearExpression
        The variable to register with the solver.
    lb : float, default=-np.inf
        Lower bound for the variable.
    ub : float, default=np.inf
        Upper bound for the variable.

    Raises
    ------
    NotImplementedError
        Must be implemented by subclasses.
    """
    raise NotImplementedError

constraints_to_solver

constraints_to_solver(constraints)

Add multiple constraints to the model.

Parameters:

Name	Type	Description	Default
constraints	list[Constraint]	List of constraints to add to the solver model.	required

Raises:

Type	Description
NotImplementedError	Must be implemented by subclasses.

Source code in fomtools/solvers/solver.py

def constraints_to_solver(self, constraints):
    """
    Add multiple constraints to the model.

    Parameters
    ----------
    constraints : list[Constraint]
        List of constraints to add to the solver model.

    Raises
    ------
    NotImplementedError
        Must be implemented by subclasses.
    """
    raise NotImplementedError

copy

copy() -> Solver

Create an independent copy of this solver.

The copy can be modified without affecting the original.

Returns:

Type	Description
Solver	A new solver instance with copied state.

Raises:

Type	Description
NotImplementedError	Must be implemented by subclasses.

Source code in fomtools/solvers/solver.py

def copy(self) -> "Solver":
    """
    Create an independent copy of this solver.

    The copy can be modified without affecting the original.

    Returns
    -------
    Solver
        A new solver instance with copied state.

    Raises
    ------
    NotImplementedError
        Must be implemented by subclasses.
    """
    raise NotImplementedError

optimize

optimize()

Solve the optimization model.

Raises:

Type	Description
NotImplementedError	Must be implemented by subclasses.

Source code in fomtools/solvers/solver.py

def optimize(self):
    """
    Solve the optimization model.

    Raises
    ------
    NotImplementedError
        Must be implemented by subclasses.
    """
    raise NotImplementedError

set_param

set_param(param, value)

Set an individual solver parameter.

Parameters:

Name	Type	Description	Default
param	str	Name of the solver parameter.	required
value	int or float	Value to assign to the parameter.	required

Raises:

Type	Description
NotImplementedError	Must be implemented by subclasses.

Source code in fomtools/solvers/solver.py

def set_param(self, param, value):
    """
    Set an individual solver parameter.

    Parameters
    ----------
    param : str
        Name of the solver parameter.
    value : int or float
        Value to assign to the parameter.

    Raises
    ------
    NotImplementedError
        Must be implemented by subclasses.
    """
    raise NotImplementedError

set_solver_params

set_solver_params(solver_params)

Apply a dictionary of solver parameters.

Parameters:

Name	Type	Description	Default
solver_params	dict	Dictionary mapping parameter names to their values. Each key-value pair is passed to set_param.	required

Source code in fomtools/solvers/solver.py

def set_solver_params(self, solver_params):
    """
    Apply a dictionary of solver parameters.

    Parameters
    ----------
    solver_params : dict
        Dictionary mapping parameter names to their values.
        Each key-value pair is passed to ``set_param``.
    """
    for k, v in solver_params.items():
        self.set_param(k, v)

verbosity_params classmethod

verbosity_params(params: dict | None) -> dict | None

The verbosity option from params, if it names this backend's option.

Used to pass only the output setting of the main solver's options on to a secondary solver of the same backend (e.g. the OBBT LP).

Source code in fomtools/solvers/solver.py

@classmethod
def verbosity_params(cls, params: dict | None) -> dict | None:
    """
    The verbosity option from ``params``, if it names this backend's option.

    Used to pass only the output setting of the main solver's options on to a
    secondary solver of the same backend (e.g. the OBBT LP).
    """
    if not params or cls.VERBOSITY_KEY is None or cls.VERBOSITY_KEY not in params:
        return None
    return {cls.VERBOSITY_KEY: params[cls.VERBOSITY_KEY]}

GurobiSolver

Bases: MIPSolver

A wrapper around the Gurobi Model that bridges LinearExpression and Constraint objects to solver-ready variables, expressions, and constraints.

Attributes:

Name	Type	Description
model	Model	The underlying Gurobi model instance.
...

Source code in fomtools/solvers/gurobi.py

class GurobiSolver(MIPSolver):
    """
    A wrapper around the Gurobi ``Model`` that bridges ``LinearExpression``
    and ``Constraint`` objects to solver-ready variables, expressions, and
    constraints.

    Attributes
    ----------
    model: gp.Model
        The underlying Gurobi model instance.

    ...

    """

    VERBOSITY_KEY = "OutputFlag"
    # Defaults applied under the user's parameters. MIPFocus 2 (prove optimality)
    # closes the bound of the verification MIPs much faster than the balanced default.
    DEFAULT_PARAMS = {"MIPFocus": 2}

    def __init__(self, solver_params=None):
        """
        Initialize an instance of the `GurobiSolver` class.

        Parameters
        ----------
        solver_params
            Gurobi parameters; they override ``DEFAULT_PARAMS``.
        """
        user = solver_params or {}
        # Gurobi parameter names are case-insensitive.
        verbosity = {k: v for k, v in user.items() if k.lower() == self.VERBOSITY_KEY.lower()}
        self.model = gp.Model(env=_env())
        # Models inherit the environment's OutputFlag=0; Gurobi's default is to log.
        verbosity = verbosity or {self.VERBOSITY_KEY: 1}
        # Apply the user's verbosity first, so a silenced model does not log the defaults.
        super().__init__({**verbosity, **self.DEFAULT_PARAMS, **user})

    def __getitem__(self, item):
        return getattr(self.model, item)

    def set_param(self, param: str, val):
        """
        Configure a solver parameter, passed verbatim to Gurobi.

        Parameters
        ----------
        param: str
            The native Gurobi parameter name to set (e.g. ``TimeLimit``,
            ``MIPGap``).
        val:
            The value to assign to the specified parameter.

        Raises
        ------
        ValueError
            If Gurobi rejects the parameter name or value.
        """
        try:
            self.model.setParam(param, val)
        except gp.GurobiError as e:
            raise ValueError(f"Invalid Gurobi solver option {param!r}={val!r}: {e}") from e

    def relax_model(self):
        """
        Build and solve the continuous relaxation of the current
        model.

        This method creates a Gurobi model in which all integrality
        constraints are relaxed to continuous domains and optimizes
        it. It then clears any cached ``relaxed_model`` property
        so that future accesses re-expose the freshly computed
        object.
        """
        self._relaxed_model = self.model.relax()
        self._relaxed_model.optimize()

        self.__dict__.pop("relaxed_model", None)

    def set_starts(self, starts: list) -> None:
        """
        Provide MIP starts, one ``{leaf: value}`` map each. Gurobi tries every start
        and keeps the best feasible one; leaves a start omits are completed by Gurobi.
        """
        if not starts:
            return
        self.model.NumStart = len(starts)
        self.model.update()
        for number, start in enumerate(starts):
            self.model.Params.StartNumber = number
            for leaf, value in start.items():
                if leaf in self.vector_var_map:
                    self.vector_var_map[leaf].Start = np.asarray(value, dtype=float).reshape(-1)
        self.model.update()

    def get_best_bound(self) -> float:
        status = self.model.Status
        usable = {
            gp.GRB.OPTIMAL,
            gp.GRB.TIME_LIMIT,
            gp.GRB.NODE_LIMIT,
            gp.GRB.ITERATION_LIMIT,
            gp.GRB.SOLUTION_LIMIT,
            gp.GRB.INTERRUPTED,
            gp.GRB.USER_OBJ_LIMIT,
            gp.GRB.WORK_LIMIT,
            gp.GRB.MEM_LIMIT,
        }
        if status not in usable:
            raise RuntimeError(f"Gurobi terminated with status {status}; no usable best bound.")
        try:
            # Continuous optimal models need not expose ObjBound.
            bound = (
                self.model.ObjVal
                if status == gp.GRB.OPTIMAL and not self.model.IsMIP
                else self.model.ObjBound
            )
        except (AttributeError, gp.GurobiError) as exc:
            raise RuntimeError(f"Gurobi status {status} has no best bound.") from exc
        if not np.isfinite(bound) or abs(bound) >= gp.GRB.INFINITY:
            raise RuntimeError(f"Gurobi status {status} has no finite best bound.")
        return float(bound)

    def has_solution(self) -> bool:
        """Whether the last solve produced a feasible solution."""
        return self.model.SolCount > 0

    def relaxed_is_optimal(self) -> bool:
        """Whether the last ``relax_model`` solve ended optimal."""
        return self._relaxed_model.Status == gp.GRB.OPTIMAL

    def add_var(self, var: LinearExpression, lb: float = -np.inf, ub: float = np.inf) -> gp.Var:
        """
        Create and register a Gurobi matrix variable corresponding to the
        given ``LinearExpression``.

        This method wraps ``model.addVar`` to create one or more decision variables
        whose shape matches that of the provided ``LinearExpression``. It also records
        the mapping between the ``LinearExpression`` object and its associated Gurobi
        variable for subsequent reference. The created variable(s) have the specified
        bounds and types.

        Parameters
        ----------
        var: LinearExpression
            LinearExpression object to be registered.

        lb: float, default=-np.inf
            Lower bound for all the elements of the matrix variable.

        ub: float, default=np.inf
            Upper bound for all the elements of the matrix variable.

        Returns
        -------
        None
        """
        for k in var.decomposition_dict.keys():
            self.vector_var_map[k] = self.model.addMVar(
                var.ndim, lb=lb, ub=ub, vtype=gp.GRB.BINARY if var.is_binary else gp.GRB.CONTINUOUS
            )
            priority = getattr(var, "branch_priority", 0)
            if var.is_binary and priority:
                self.vector_var_map[k].BranchPriority = priority
        self.model.update()

    def add_param_stack_var(
        self, param_stack_var: LinearExpression, param_list: list[LinearExpression]
    ) -> None:
        """
        Create a stacked Gurobi ``MVar`` by horizontally concatenating
        existing parameter variables.

        Registers ``param_stack_var`` in ``self.vector_var_map`` as the
        horizontal stack of the variables associated with each element of
        ``param_list``. This is useful for building block-structured vectors
        from previously created components.

        Parameters
        ----------
        param_stack_var: LinearExpression
            LinearExpression object to be registered.

        param_list: list[LinearExpression]
            List of LinearExpression objects whose associated Gurobi variables are to be stacked.
            Each entry must already exist in ``self.vector_var_map`` and map to a compatible
            ``gp.Var`` or ``gp.MVar`` object.

        Returns
        -------
        None

        Raises
        ------
        KeyError
            If any element of ``param_list`` is not already registered to the model.

        """
        try:
            self.vector_var_map[param_stack_var] = gp.hstack(
                [self.vector_var_map[param] for param in param_list]
            )
        except KeyError:
            raise KeyError(
                """One or more elements in 'param_list' are not registered with the Gurobi model.
                Ensure that all parameters in 'param_list' have been added to 'vector_var_map' before stacking."""
            )

    def leaf_expr_to_solver(self, expression: LinearExpression) -> LinearExpression:
        """
        Resolve a leaf ``LinearExpression`` into its solver-side Gurobi variable.

        Parameters
        ----------
        expression: LinearExpression
            A leaf ``LinearExpression`` object with a single entry in
            ``expression.decomposition_dict``.

        Returns
        -------
        gp.Var or gp.MVar
            The Gurobi variable corresponding to the leaf ``expression``'s sole key.

        Raises
        ------
        TypeError
            If ``expression`` is not a leaf ``LinearExpression``.
        """
        if expression.is_leaf:
            return self.vector_var_map[next(iter(expression.decomposition_dict))]

        raise TypeError(
            "Invalid input: expression must be a leaf. Received a composite one instead."
        )

    def expression_to_solver(self, expression: LinearExpression) -> gp.MLinExpr:
        r"""
        Convert a ``LinearExpression`` into a Gurobi ``MLinExpr``.

        Given a ``LinearExpression`` object with a decomposition of the form
        ``expression.decomposition_dict = {key_1: matrix_1, ..., key_n: matrix_n}``,
        this method constructs the equivalent solver-side representation by multiplying each
        coefficient ``matrix_i`` with its corresponding Gurobi variable and
        summing the results:

        $$
            \sum_{i}^{N} matrix_i @ variable_i
        $$

        The resulting object is a Gurobi linear expression with the same dimensionality as
        the input ``expression``, and can be directly used in model constraints or
        objective functions.

        Parameters
        ----------
        expression: LinearExpression
            The LinearExpression object to be converted into its corresponding Gurobi representation.

        Returns
        -------
        gp.MLinExpr
            The Gurobi linear expression corresponding to the input ``expression``, representing
            its solver-side linear form.
        """
        gp_expr = np.zeros(expression.ndim)
        for k, v in expression.decomposition_dict.items():
            if self.vector_var_map.get(k, None) is None:
                self.add_var(k)
            gp_expr += v @ self.vector_var_map[k]

        return gp_expr + expression.const

    def constraints_to_solver(self, constraints: List[Constraint]):
        """
        Add multiple constraints to the Gurobi model.

        Parameters
        ----------
        constraints: list[Constraint]
            A list of constraint expressions, typically of the form ``lhs == rhs``,
            ``lhs <= rhs``, or ``lhs >= rhs`` where ``lhs`` and ``rhs`` are of type
            ``LinearExpression``.

        Returns
        -------
        None

        """
        for constr in constraints:
            self.add_constr(constr)

    def update_bounds(
        self, expression: LinearExpression, lb: float | np.ndarray, ub: float | np.ndarray
    ):
        """
        Update the lower and upper bounds of the solver variables corresponding
        to the given ``LinearExpression``.

        Parameters
        ----------
        expression: LinearExpression
            A ``LinearExpression`` object registered with the solver whose
            underlying variable bounds are to be updated.

        lb: float or np.ndarray
            New lower bound(s).

        ub: float or np.ndarray
            New upper bound(s).

        """
        var = self.expression_to_solver(expression)
        var.LB = lb
        var.UB = ub

    def add_constr(self, constraint) -> None:
        """
        Add a single constraint to the Gurobi model.

        Parameters
        ----------
        constraint: Constraint | gp.TempConstr
            A compatible Gurobi constraint expression. Typical forms include
            comparisons like ``lhs == rhs``, ``lhs <= rhs``, where ``lhs``,
            and ``rhs`` are of type ``LinearExpression`` or ``gp.MLinExpr``.

        Returns
        -------
        None
        """
        try:
            self.model.addConstr(constraint)
        except gp.GurobiError:
            # Expression on the right-hand side of the constraint object is always 0.
            lhs = self.expression_to_solver(constraint.lhs)
            match constraint.operator:
                case "=":
                    self.model.addConstr(lhs == 0)
                case "<=":
                    self.model.addConstr(lhs <= 0)
            self.model.update()

    def set_objective(self, expression: LinearExpression, sense: str = None) -> None:
        """
        Set the optimization objective for the underlying Gurobi model.

        Parameters
        ----------
        expression: LinearExpression or 0
            A ``LinearExpression`` object registered with the solver representing the
            optimization objective, or 0 if setting a zero objective.

        sense: {"minimize", "maximize"}
            Objective sense specified as a string.

        Returns
        -------
        None

        Raises
        ------
        TypeError
            If provided objective ``sense`` is not ``"minimize"`` or ``"maximize"``.

        ValueError
            If provided expression is anything other than a ``LinearExpression`` object or a 0.

        """
        match sense:
            case "minimize":
                sense = gp.GRB.MINIMIZE
            case "maximize":
                sense = gp.GRB.MAXIMIZE
            case None:
                sense = None
            case _:
                raise TypeError(f"Objective sense must be either minimize or maximize, got {sense}")

        try:
            if int(expression) != 0:
                raise ValueError(
                    f"Objective function must be either a LinearExpression object or a 0, got {expression}"
                )

            self.model.setObjective(expression, sense=sense)
        except TypeError:
            self.model.setObjective(self.expression_to_solver(expression), sense=sense)
        self.model.update()

    def add_sum_squares_constr(self, expr, radius_sq, sense) -> None:
        """
        Add a sum-of-squares (quadratic) constraint to the Gurobi model.

        Encodes ``sum_i (a_i^T x + b_i)^2  op  radius_sq`` where ``op`` is given by
        ``sense``. The squared two-norm is built as ``gp_expr @ gp_expr`` (a dot product),
        the same construct used by `set_two_norm_objective`. The const term ``b`` is
        included via `expression_to_solver`.

        The ``>=`` / ``==`` senses are nonconvex quadratic constraints; ``NonConvex=2`` is
        set for them so Gurobi accepts the model across versions (the convex ``<=`` sense
        does not need it).

        Parameters
        ----------
        expr : LinearExpression
            The expression whose squared two-norm is constrained.
        radius_sq : float
            The nonnegative scalar right-hand side (e.g. ``r**2``).
        sense : str {"<=", ">=", "=="}
            The relational operator.
        """
        gp_expr = self.expression_to_solver(expr)
        ss = gp_expr @ gp_expr  # 0-d MQuadExpr = sum of squares (matrix API)
        if sense != "<=":
            self.model.Params.NonConvex = 2
        # The matrix-API MQuadExpr is added via addConstr (addQConstr rejects the 0-d MQuadExpr).
        match sense:
            case "<=":
                self.model.addConstr(ss <= radius_sq)
            case ">=":
                self.model.addConstr(ss >= radius_sq)
            case "==":
                self.model.addConstr(ss == radius_sq)
        self.model.update()

    def set_two_norm_objective(self, expr_list, lb_list, ub_list) -> None:
        """
        Set a quadratic two-norm objective for the underlying Gurobi model.

        The objective is to maximize sum_i ||expr_i||_2^2 = sum_i (expr_i @ expr_i).

        Parameters
        ----------
        expr_list : list[LinearExpression]
            List of expressions whose squared norms are summed in the objective.
        lb_list : list[np.ndarray]
            Lower bounds for each expression (currently unused but passed for interface consistency).
        ub_list : list[np.ndarray]
            Upper bounds for each expression (currently unused but passed for interface consistency).

        Returns
        -------
        None
        """
        quadratic_obj = 0
        for expr in expr_list:
            gp_expr = self.expression_to_solver(expr)
            quadratic_obj = quadratic_obj + gp_expr @ gp_expr

        self.model.setObjective(quadratic_obj, sense=gp.GRB.MAXIMIZE)
        self.model.update()

    def optimize(self):
        """
        Optimize the underlying Gurobi model.

        Notes
        -----
        The algorithm used depends on the model objective and
        on the parameter settings
        """
        self.model.optimize()

    def get_solution_value(self, expression: LinearExpression) -> np.ndarray:
        """
        Compute the solution value for the given ``LinearExpression``
        object.

        For each leaf expression in the decomposition of the given
        ``LinearExpression`` object, this method retrieves the corresponding
        solver variables from the model and computes the dot product between
        the coefficient matrix and the vector of variable values associated with
        that leaf.

        Parameters
        ----------
        expression: LinearExpression
            A ``LinearExpression`` object registered with the solver. Each element in
            its decomposition must have a corresponding entry in ``self.vector_var_map``.

        Returns
        -------
        np.ndarray[float]
            A 1-D array of shape ``(expression.ndim,)`` containing the evaluated values
            of the expression at the current model solution.
        """
        if not self.has_solution():
            raise RuntimeError("No feasible incumbent available.")
        out = np.array(expression.const, dtype=float, copy=True)
        for k, v in expression.decomposition_dict.items():
            out += v @ [self.model.getVarByName(var).X for var in self.vector_var_map[k].VarName]
        return out

    def get_relaxed_variables(self, expression: LinearExpression) -> Tuple[gp.MVar, np.ndarray]:
        """
        Reconstruct solver-side variables from a ``LinearExpression`` object
        using the relaxed Gurobi model.

        Parameters
        ----------
        expression: LinearExpression
            A ``LinearExpression`` object registered with the solver. Each element in
            its decomposition must have a corresponding entry in ``self.vector_var_map``.

        Returns
        -------
        Tuple[gp.MVar, np.ndarray]
            A single Gurobi matrix variable formed by horizontally stacking the
            individual ``gp.MVar`` objects corresponding to each element of
            the input ``LinearExpression``, and a NumPy array containing the
            current values of those stacked variables.
        """
        # If no relaxed model currently exists, build and solve it first.
        if self.relaxed_model is None:
            self.relax_model()

        gp_vars = []
        for k in expression.decomposition_dict.keys():
            gp_vars.append(
                gp.MVar.fromlist(
                    [self.relaxed_model.getVarByName(i) for i in self.vector_var_map[k].VarName]
                )
            )
        return gp.hstack(gp_vars), gp.hstack(gp_vars).X

    def is_optimal(self) -> bool:
        """
        Check whether the underlying Gurobi model has reached an optimal solution.
        """
        return self["status"] == gp.GRB.OPTIMAL

    def copy(self) -> "GurobiSolver":
        """
        Create an independent copy of this Gurobi solver.

        Returns a new GurobiSolver with a copied model and rebuilt
        vector_var_map pointing to variables in the copied model.

        Returns
        -------
        GurobiSolver
            A new solver instance with copied state.
        """
        new_solver = GurobiSolver.__new__(GurobiSolver)
        new_solver._relaxed_model = None
        new_solver.solver_params = self.solver_params
        new_solver.model = self.model.copy()

        # Rebuild vector_var_map with variables from the copied model
        new_solver.vector_var_map = {}
        for k, v in self.vector_var_map.items():
            # v is a gp.MVar - reconstruct it from variable names
            new_solver.vector_var_map[k] = gp.MVar.fromlist(
                [new_solver.model.getVarByName(name) for name in v.VarName]
            )

        return new_solver

__init__

__init__(solver_params=None)

Parameters:

Name	Type	Description	Default
solver_params		Gurobi parameters; they override DEFAULT_PARAMS.	None

Source code in fomtools/solvers/gurobi.py

def __init__(self, solver_params=None):
    """
    Initialize an instance of the `GurobiSolver` class.

    Parameters
    ----------
    solver_params
        Gurobi parameters; they override ``DEFAULT_PARAMS``.
    """
    user = solver_params or {}
    # Gurobi parameter names are case-insensitive.
    verbosity = {k: v for k, v in user.items() if k.lower() == self.VERBOSITY_KEY.lower()}
    self.model = gp.Model(env=_env())
    # Models inherit the environment's OutputFlag=0; Gurobi's default is to log.
    verbosity = verbosity or {self.VERBOSITY_KEY: 1}
    # Apply the user's verbosity first, so a silenced model does not log the defaults.
    super().__init__({**verbosity, **self.DEFAULT_PARAMS, **user})

add_constr

add_constr(constraint) -> None

Add a single constraint to the Gurobi model.

Parameters:

Name	Type	Description	Default
constraint		A compatible Gurobi constraint expression. Typical forms include comparisons like lhs == rhs, lhs <= rhs, where lhs, and rhs are of type LinearExpression or gp.MLinExpr.	required

Returns:

Type	Description
None

Source code in fomtools/solvers/gurobi.py

def add_constr(self, constraint) -> None:
    """
    Add a single constraint to the Gurobi model.

    Parameters
    ----------
    constraint: Constraint | gp.TempConstr
        A compatible Gurobi constraint expression. Typical forms include
        comparisons like ``lhs == rhs``, ``lhs <= rhs``, where ``lhs``,
        and ``rhs`` are of type ``LinearExpression`` or ``gp.MLinExpr``.

    Returns
    -------
    None
    """
    try:
        self.model.addConstr(constraint)
    except gp.GurobiError:
        # Expression on the right-hand side of the constraint object is always 0.
        lhs = self.expression_to_solver(constraint.lhs)
        match constraint.operator:
            case "=":
                self.model.addConstr(lhs == 0)
            case "<=":
                self.model.addConstr(lhs <= 0)
        self.model.update()

add_param_stack_var

add_param_stack_var(param_stack_var: LinearExpression, param_list: list[LinearExpression]) -> None

Create a stacked Gurobi MVar by horizontally concatenating existing parameter variables.

Registers param_stack_var in self.vector_var_map as the horizontal stack of the variables associated with each element of param_list. This is useful for building block-structured vectors from previously created components.

Parameters:

Name	Type	Description	Default
param_stack_var	LinearExpression	LinearExpression object to be registered.	required
param_list	list[LinearExpression]	List of LinearExpression objects whose associated Gurobi variables are to be stacked. Each entry must already exist in self.vector_var_map and map to a compatible gp.Var or gp.MVar object.	required

Returns:

Type	Description
None

Raises:

Type	Description
KeyError	If any element of param_list is not already registered to the model.

Source code in fomtools/solvers/gurobi.py

def add_param_stack_var(
    self, param_stack_var: LinearExpression, param_list: list[LinearExpression]
) -> None:
    """
    Create a stacked Gurobi ``MVar`` by horizontally concatenating
    existing parameter variables.

    Registers ``param_stack_var`` in ``self.vector_var_map`` as the
    horizontal stack of the variables associated with each element of
    ``param_list``. This is useful for building block-structured vectors
    from previously created components.

    Parameters
    ----------
    param_stack_var: LinearExpression
        LinearExpression object to be registered.

    param_list: list[LinearExpression]
        List of LinearExpression objects whose associated Gurobi variables are to be stacked.
        Each entry must already exist in ``self.vector_var_map`` and map to a compatible
        ``gp.Var`` or ``gp.MVar`` object.

    Returns
    -------
    None

    Raises
    ------
    KeyError
        If any element of ``param_list`` is not already registered to the model.

    """
    try:
        self.vector_var_map[param_stack_var] = gp.hstack(
            [self.vector_var_map[param] for param in param_list]
        )
    except KeyError:
        raise KeyError(
            """One or more elements in 'param_list' are not registered with the Gurobi model.
            Ensure that all parameters in 'param_list' have been added to 'vector_var_map' before stacking."""
        )

add_sum_squares_constr

add_sum_squares_constr(expr, radius_sq, sense) -> None

Add a sum-of-squares (quadratic) constraint to the Gurobi model.

Encodes sum_i (a_i^T x + b_i)^2 op radius_sq where op is given by sense. The squared two-norm is built as gp_expr @ gp_expr (a dot product), the same construct used by set_two_norm_objective. The const term b is included via expression_to_solver.

The >= / == senses are nonconvex quadratic constraints; NonConvex=2 is set for them so Gurobi accepts the model across versions (the convex <= sense does not need it).

Parameters:

Name	Type	Description	Default
expr	LinearExpression	The expression whose squared two-norm is constrained.	required
radius_sq	float	The nonnegative scalar right-hand side (e.g. r**2).	required
sense	str {"<=", ">=", "=="}	The relational operator.	required

Source code in fomtools/solvers/gurobi.py

def add_sum_squares_constr(self, expr, radius_sq, sense) -> None:
    """
    Add a sum-of-squares (quadratic) constraint to the Gurobi model.

    Encodes ``sum_i (a_i^T x + b_i)^2  op  radius_sq`` where ``op`` is given by
    ``sense``. The squared two-norm is built as ``gp_expr @ gp_expr`` (a dot product),
    the same construct used by `set_two_norm_objective`. The const term ``b`` is
    included via `expression_to_solver`.

    The ``>=`` / ``==`` senses are nonconvex quadratic constraints; ``NonConvex=2`` is
    set for them so Gurobi accepts the model across versions (the convex ``<=`` sense
    does not need it).

    Parameters
    ----------
    expr : LinearExpression
        The expression whose squared two-norm is constrained.
    radius_sq : float
        The nonnegative scalar right-hand side (e.g. ``r**2``).
    sense : str {"<=", ">=", "=="}
        The relational operator.
    """
    gp_expr = self.expression_to_solver(expr)
    ss = gp_expr @ gp_expr  # 0-d MQuadExpr = sum of squares (matrix API)
    if sense != "<=":
        self.model.Params.NonConvex = 2
    # The matrix-API MQuadExpr is added via addConstr (addQConstr rejects the 0-d MQuadExpr).
    match sense:
        case "<=":
            self.model.addConstr(ss <= radius_sq)
        case ">=":
            self.model.addConstr(ss >= radius_sq)
        case "==":
            self.model.addConstr(ss == radius_sq)
    self.model.update()

add_var

add_var(var: LinearExpression, lb: float = -inf, ub: float = inf) -> Var

Create and register a Gurobi matrix variable corresponding to the given LinearExpression.

This method wraps model.addVar to create one or more decision variables whose shape matches that of the provided LinearExpression. It also records the mapping between the LinearExpression object and its associated Gurobi variable for subsequent reference. The created variable(s) have the specified bounds and types.

Parameters:

Name	Type	Description	Default
var	LinearExpression	LinearExpression object to be registered.	required
lb	float	Lower bound for all the elements of the matrix variable.	-inf
ub	float	Upper bound for all the elements of the matrix variable.	inf

Returns:

Type	Description
None

Source code in fomtools/solvers/gurobi.py

def add_var(self, var: LinearExpression, lb: float = -np.inf, ub: float = np.inf) -> gp.Var:
    """
    Create and register a Gurobi matrix variable corresponding to the
    given ``LinearExpression``.

    This method wraps ``model.addVar`` to create one or more decision variables
    whose shape matches that of the provided ``LinearExpression``. It also records
    the mapping between the ``LinearExpression`` object and its associated Gurobi
    variable for subsequent reference. The created variable(s) have the specified
    bounds and types.

    Parameters
    ----------
    var: LinearExpression
        LinearExpression object to be registered.

    lb: float, default=-np.inf
        Lower bound for all the elements of the matrix variable.

    ub: float, default=np.inf
        Upper bound for all the elements of the matrix variable.

    Returns
    -------
    None
    """
    for k in var.decomposition_dict.keys():
        self.vector_var_map[k] = self.model.addMVar(
            var.ndim, lb=lb, ub=ub, vtype=gp.GRB.BINARY if var.is_binary else gp.GRB.CONTINUOUS
        )
        priority = getattr(var, "branch_priority", 0)
        if var.is_binary and priority:
            self.vector_var_map[k].BranchPriority = priority
    self.model.update()

constraints_to_solver

constraints_to_solver(constraints: List[Constraint])

Add multiple constraints to the Gurobi model.

Parameters:

Name	Type	Description	Default
constraints	List[Constraint]	A list of constraint expressions, typically of the form lhs == rhs, lhs <= rhs, or lhs >= rhs where lhs and rhs are of type LinearExpression.	required

Returns:

Type	Description
None

Source code in fomtools/solvers/gurobi.py

def constraints_to_solver(self, constraints: List[Constraint]):
    """
    Add multiple constraints to the Gurobi model.

    Parameters
    ----------
    constraints: list[Constraint]
        A list of constraint expressions, typically of the form ``lhs == rhs``,
        ``lhs <= rhs``, or ``lhs >= rhs`` where ``lhs`` and ``rhs`` are of type
        ``LinearExpression``.

    Returns
    -------
    None

    """
    for constr in constraints:
        self.add_constr(constr)

copy

copy() -> GurobiSolver

Create an independent copy of this Gurobi solver.

Returns a new GurobiSolver with a copied model and rebuilt vector_var_map pointing to variables in the copied model.

Returns:

Type	Description
GurobiSolver	A new solver instance with copied state.

Source code in fomtools/solvers/gurobi.py

def copy(self) -> "GurobiSolver":
    """
    Create an independent copy of this Gurobi solver.

    Returns a new GurobiSolver with a copied model and rebuilt
    vector_var_map pointing to variables in the copied model.

    Returns
    -------
    GurobiSolver
        A new solver instance with copied state.
    """
    new_solver = GurobiSolver.__new__(GurobiSolver)
    new_solver._relaxed_model = None
    new_solver.solver_params = self.solver_params
    new_solver.model = self.model.copy()

    # Rebuild vector_var_map with variables from the copied model
    new_solver.vector_var_map = {}
    for k, v in self.vector_var_map.items():
        # v is a gp.MVar - reconstruct it from variable names
        new_solver.vector_var_map[k] = gp.MVar.fromlist(
            [new_solver.model.getVarByName(name) for name in v.VarName]
        )

    return new_solver

expression_to_solver

expression_to_solver(expression: LinearExpression) -> MLinExpr

Convert a LinearExpression into a Gurobi MLinExpr.

Given a LinearExpression object with a decomposition of the form expression.decomposition_dict = {key_1: matrix_1, ..., key_n: matrix_n}, this method constructs the equivalent solver-side representation by multiplying each coefficient matrix_i with its corresponding Gurobi variable and summing the results:

\[ \sum_{i}^{N} matrix_i @ variable_i \]

The resulting object is a Gurobi linear expression with the same dimensionality as the input expression, and can be directly used in model constraints or objective functions.

Parameters:

Name	Type	Description	Default
expression	LinearExpression	The LinearExpression object to be converted into its corresponding Gurobi representation.	required

Returns:

Type	Description
MLinExpr	The Gurobi linear expression corresponding to the input expression, representing its solver-side linear form.

Source code in fomtools/solvers/gurobi.py

def expression_to_solver(self, expression: LinearExpression) -> gp.MLinExpr:
    r"""
    Convert a ``LinearExpression`` into a Gurobi ``MLinExpr``.

    Given a ``LinearExpression`` object with a decomposition of the form
    ``expression.decomposition_dict = {key_1: matrix_1, ..., key_n: matrix_n}``,
    this method constructs the equivalent solver-side representation by multiplying each
    coefficient ``matrix_i`` with its corresponding Gurobi variable and
    summing the results:

    $$
        \sum_{i}^{N} matrix_i @ variable_i
    $$

    The resulting object is a Gurobi linear expression with the same dimensionality as
    the input ``expression``, and can be directly used in model constraints or
    objective functions.

    Parameters
    ----------
    expression: LinearExpression
        The LinearExpression object to be converted into its corresponding Gurobi representation.

    Returns
    -------
    gp.MLinExpr
        The Gurobi linear expression corresponding to the input ``expression``, representing
        its solver-side linear form.
    """
    gp_expr = np.zeros(expression.ndim)
    for k, v in expression.decomposition_dict.items():
        if self.vector_var_map.get(k, None) is None:
            self.add_var(k)
        gp_expr += v @ self.vector_var_map[k]

    return gp_expr + expression.const

get_relaxed_variables

get_relaxed_variables(expression: LinearExpression) -> Tuple[MVar, ndarray]

Reconstruct solver-side variables from a LinearExpression object using the relaxed Gurobi model.

Parameters:

Name	Type	Description	Default
expression	LinearExpression	A LinearExpression object registered with the solver. Each element in its decomposition must have a corresponding entry in self.vector_var_map.	required

Returns:

Type	Description
Tuple[MVar, ndarray]	A single Gurobi matrix variable formed by horizontally stacking the individual gp.MVar objects corresponding to each element of the input LinearExpression, and a NumPy array containing the current values of those stacked variables.

Source code in fomtools/solvers/gurobi.py

def get_relaxed_variables(self, expression: LinearExpression) -> Tuple[gp.MVar, np.ndarray]:
    """
    Reconstruct solver-side variables from a ``LinearExpression`` object
    using the relaxed Gurobi model.

    Parameters
    ----------
    expression: LinearExpression
        A ``LinearExpression`` object registered with the solver. Each element in
        its decomposition must have a corresponding entry in ``self.vector_var_map``.

    Returns
    -------
    Tuple[gp.MVar, np.ndarray]
        A single Gurobi matrix variable formed by horizontally stacking the
        individual ``gp.MVar`` objects corresponding to each element of
        the input ``LinearExpression``, and a NumPy array containing the
        current values of those stacked variables.
    """
    # If no relaxed model currently exists, build and solve it first.
    if self.relaxed_model is None:
        self.relax_model()

    gp_vars = []
    for k in expression.decomposition_dict.keys():
        gp_vars.append(
            gp.MVar.fromlist(
                [self.relaxed_model.getVarByName(i) for i in self.vector_var_map[k].VarName]
            )
        )
    return gp.hstack(gp_vars), gp.hstack(gp_vars).X

get_solution_value

get_solution_value(expression: LinearExpression) -> ndarray

Compute the solution value for the given LinearExpression object.

For each leaf expression in the decomposition of the given LinearExpression object, this method retrieves the corresponding solver variables from the model and computes the dot product between the coefficient matrix and the vector of variable values associated with that leaf.

Parameters:

Name	Type	Description	Default
expression	LinearExpression	A LinearExpression object registered with the solver. Each element in its decomposition must have a corresponding entry in self.vector_var_map.	required

Returns:

Type	Description
ndarray[float]	A 1-D array of shape (expression.ndim,) containing the evaluated values of the expression at the current model solution.

Source code in fomtools/solvers/gurobi.py

def get_solution_value(self, expression: LinearExpression) -> np.ndarray:
    """
    Compute the solution value for the given ``LinearExpression``
    object.

    For each leaf expression in the decomposition of the given
    ``LinearExpression`` object, this method retrieves the corresponding
    solver variables from the model and computes the dot product between
    the coefficient matrix and the vector of variable values associated with
    that leaf.

    Parameters
    ----------
    expression: LinearExpression
        A ``LinearExpression`` object registered with the solver. Each element in
        its decomposition must have a corresponding entry in ``self.vector_var_map``.

    Returns
    -------
    np.ndarray[float]
        A 1-D array of shape ``(expression.ndim,)`` containing the evaluated values
        of the expression at the current model solution.
    """
    if not self.has_solution():
        raise RuntimeError("No feasible incumbent available.")
    out = np.array(expression.const, dtype=float, copy=True)
    for k, v in expression.decomposition_dict.items():
        out += v @ [self.model.getVarByName(var).X for var in self.vector_var_map[k].VarName]
    return out

has_solution

has_solution() -> bool

Whether the last solve produced a feasible solution.

Source code in fomtools/solvers/gurobi.py

def has_solution(self) -> bool:
    """Whether the last solve produced a feasible solution."""
    return self.model.SolCount > 0

is_optimal

is_optimal() -> bool

Check whether the underlying Gurobi model has reached an optimal solution.

Source code in fomtools/solvers/gurobi.py

def is_optimal(self) -> bool:
    """
    Check whether the underlying Gurobi model has reached an optimal solution.
    """
    return self["status"] == gp.GRB.OPTIMAL

leaf_expr_to_solver

leaf_expr_to_solver(expression: LinearExpression) -> LinearExpression

Resolve a leaf LinearExpression into its solver-side Gurobi variable.

Parameters:

Name	Type	Description	Default
expression	LinearExpression	A leaf LinearExpression object with a single entry in expression.decomposition_dict.	required

Returns:

Type	Description
Var or MVar	The Gurobi variable corresponding to the leaf expression's sole key.

Raises:

Type	Description
TypeError	If expression is not a leaf LinearExpression.

Source code in fomtools/solvers/gurobi.py

def leaf_expr_to_solver(self, expression: LinearExpression) -> LinearExpression:
    """
    Resolve a leaf ``LinearExpression`` into its solver-side Gurobi variable.

    Parameters
    ----------
    expression: LinearExpression
        A leaf ``LinearExpression`` object with a single entry in
        ``expression.decomposition_dict``.

    Returns
    -------
    gp.Var or gp.MVar
        The Gurobi variable corresponding to the leaf ``expression``'s sole key.

    Raises
    ------
    TypeError
        If ``expression`` is not a leaf ``LinearExpression``.
    """
    if expression.is_leaf:
        return self.vector_var_map[next(iter(expression.decomposition_dict))]

    raise TypeError(
        "Invalid input: expression must be a leaf. Received a composite one instead."
    )

optimize

optimize()

Optimize the underlying Gurobi model.

Notes

The algorithm used depends on the model objective and on the parameter settings

Source code in fomtools/solvers/gurobi.py

def optimize(self):
    """
    Optimize the underlying Gurobi model.

    Notes
    -----
    The algorithm used depends on the model objective and
    on the parameter settings
    """
    self.model.optimize()

relax_model

relax_model()

Build and solve the continuous relaxation of the current model.

This method creates a Gurobi model in which all integrality constraints are relaxed to continuous domains and optimizes it. It then clears any cached relaxed_model property so that future accesses re-expose the freshly computed object.

Source code in fomtools/solvers/gurobi.py

def relax_model(self):
    """
    Build and solve the continuous relaxation of the current
    model.

    This method creates a Gurobi model in which all integrality
    constraints are relaxed to continuous domains and optimizes
    it. It then clears any cached ``relaxed_model`` property
    so that future accesses re-expose the freshly computed
    object.
    """
    self._relaxed_model = self.model.relax()
    self._relaxed_model.optimize()

    self.__dict__.pop("relaxed_model", None)

relaxed_is_optimal

relaxed_is_optimal() -> bool

Whether the last relax_model solve ended optimal.

Source code in fomtools/solvers/gurobi.py

def relaxed_is_optimal(self) -> bool:
    """Whether the last ``relax_model`` solve ended optimal."""
    return self._relaxed_model.Status == gp.GRB.OPTIMAL

set_objective

set_objective(expression: LinearExpression, sense: str = None) -> None

Set the optimization objective for the underlying Gurobi model.

Parameters:

Name	Type	Description	Default
expression	LinearExpression	A LinearExpression object registered with the solver representing the optimization objective, or 0 if setting a zero objective.	required
sense	str	Objective sense specified as a string.	None

Returns:

Type	Description
None

Raises:

Type	Description
TypeError	If provided objective sense is not "minimize" or "maximize".
ValueError	If provided expression is anything other than a LinearExpression object or a 0.

Source code in fomtools/solvers/gurobi.py

def set_objective(self, expression: LinearExpression, sense: str = None) -> None:
    """
    Set the optimization objective for the underlying Gurobi model.

    Parameters
    ----------
    expression: LinearExpression or 0
        A ``LinearExpression`` object registered with the solver representing the
        optimization objective, or 0 if setting a zero objective.

    sense: {"minimize", "maximize"}
        Objective sense specified as a string.

    Returns
    -------
    None

    Raises
    ------
    TypeError
        If provided objective ``sense`` is not ``"minimize"`` or ``"maximize"``.

    ValueError
        If provided expression is anything other than a ``LinearExpression`` object or a 0.

    """
    match sense:
        case "minimize":
            sense = gp.GRB.MINIMIZE
        case "maximize":
            sense = gp.GRB.MAXIMIZE
        case None:
            sense = None
        case _:
            raise TypeError(f"Objective sense must be either minimize or maximize, got {sense}")

    try:
        if int(expression) != 0:
            raise ValueError(
                f"Objective function must be either a LinearExpression object or a 0, got {expression}"
            )

        self.model.setObjective(expression, sense=sense)
    except TypeError:
        self.model.setObjective(self.expression_to_solver(expression), sense=sense)
    self.model.update()

set_param

set_param(param: str, val)

Configure a solver parameter, passed verbatim to Gurobi.

Parameters:

Name	Type	Description	Default
param	str	The native Gurobi parameter name to set (e.g. TimeLimit, MIPGap).	required
val		The value to assign to the specified parameter.	required

Raises:

Type	Description
ValueError	If Gurobi rejects the parameter name or value.

Source code in fomtools/solvers/gurobi.py

def set_param(self, param: str, val):
    """
    Configure a solver parameter, passed verbatim to Gurobi.

    Parameters
    ----------
    param: str
        The native Gurobi parameter name to set (e.g. ``TimeLimit``,
        ``MIPGap``).
    val:
        The value to assign to the specified parameter.

    Raises
    ------
    ValueError
        If Gurobi rejects the parameter name or value.
    """
    try:
        self.model.setParam(param, val)
    except gp.GurobiError as e:
        raise ValueError(f"Invalid Gurobi solver option {param!r}={val!r}: {e}") from e

set_starts

set_starts(starts: list) -> None

Provide MIP starts, one {leaf: value} map each. Gurobi tries every start and keeps the best feasible one; leaves a start omits are completed by Gurobi.

Source code in fomtools/solvers/gurobi.py

def set_starts(self, starts: list) -> None:
    """
    Provide MIP starts, one ``{leaf: value}`` map each. Gurobi tries every start
    and keeps the best feasible one; leaves a start omits are completed by Gurobi.
    """
    if not starts:
        return
    self.model.NumStart = len(starts)
    self.model.update()
    for number, start in enumerate(starts):
        self.model.Params.StartNumber = number
        for leaf, value in start.items():
            if leaf in self.vector_var_map:
                self.vector_var_map[leaf].Start = np.asarray(value, dtype=float).reshape(-1)
    self.model.update()

set_two_norm_objective

set_two_norm_objective(expr_list, lb_list, ub_list) -> None

Set a quadratic two-norm objective for the underlying Gurobi model.

The objective is to maximize sum_i ||expr_i||_2^2 = sum_i (expr_i @ expr_i).

Parameters:

Name	Type	Description	Default
expr_list	list[LinearExpression]	List of expressions whose squared norms are summed in the objective.	required
lb_list	list[ndarray]	Lower bounds for each expression (currently unused but passed for interface consistency).	required
ub_list	list[ndarray]	Upper bounds for each expression (currently unused but passed for interface consistency).	required

Returns:

Type	Description
None

Source code in fomtools/solvers/gurobi.py

def set_two_norm_objective(self, expr_list, lb_list, ub_list) -> None:
    """
    Set a quadratic two-norm objective for the underlying Gurobi model.

    The objective is to maximize sum_i ||expr_i||_2^2 = sum_i (expr_i @ expr_i).

    Parameters
    ----------
    expr_list : list[LinearExpression]
        List of expressions whose squared norms are summed in the objective.
    lb_list : list[np.ndarray]
        Lower bounds for each expression (currently unused but passed for interface consistency).
    ub_list : list[np.ndarray]
        Upper bounds for each expression (currently unused but passed for interface consistency).

    Returns
    -------
    None
    """
    quadratic_obj = 0
    for expr in expr_list:
        gp_expr = self.expression_to_solver(expr)
        quadratic_obj = quadratic_obj + gp_expr @ gp_expr

    self.model.setObjective(quadratic_obj, sense=gp.GRB.MAXIMIZE)
    self.model.update()

update_bounds

update_bounds(expression: LinearExpression, lb: float | ndarray, ub: float | ndarray)

Update the lower and upper bounds of the solver variables corresponding to the given LinearExpression.

Parameters:

Name	Type	Description	Default
expression	LinearExpression	A LinearExpression object registered with the solver whose underlying variable bounds are to be updated.	required
lb	float | ndarray	New lower bound(s).	required
ub	float | ndarray	New upper bound(s).	required

Source code in fomtools/solvers/gurobi.py

def update_bounds(
    self, expression: LinearExpression, lb: float | np.ndarray, ub: float | np.ndarray
):
    """
    Update the lower and upper bounds of the solver variables corresponding
    to the given ``LinearExpression``.

    Parameters
    ----------
    expression: LinearExpression
        A ``LinearExpression`` object registered with the solver whose
        underlying variable bounds are to be updated.

    lb: float or np.ndarray
        New lower bound(s).

    ub: float or np.ndarray
        New upper bound(s).

    """
    var = self.expression_to_solver(expression)
    var.LB = lb
    var.UB = ub

HiGHSSolver

Bases: MIPSolver

Source code in fomtools/solvers/highs.py

class HiGHSSolver(MIPSolver):
    VERBOSITY_KEY = "output_flag"

    def __init__(self, solver_params=None):
        # Variable to keep track of the number of HiGHS variables created.
        self._var_counter = -1
        self.model = hp.Highs()
        super().__init__(solver_params)

    def __getitem__(self, item):
        if item == "objBound":
            return self.get_best_bound()
        if item == "status":
            return self.model.getModelStatus()
        if item == "objVal":
            if not self.has_solution():
                raise RuntimeError("HiGHS has no feasible incumbent.")
            return self.model.getInfo().objective_function_value
        raise KeyError(item)

    def set_param(self, param: str, val):
        """
        Configure a solver option, passed verbatim to HiGHS.

        The option name must be a native HiGHS option (e.g. ``time_limit``,
        ``mip_rel_gap``, ``output_flag``); no name translation is performed.

        Parameters
        ----------
        param: str
            The native HiGHS option name to set.
        val:
            The value to assign to the specified option.

        Raises
        ------
        ValueError
            If HiGHS rejects the option name or value.
        """
        status = self.model.setOptionValue(param, val)
        if status == hp.HighsStatus.kError:
            raise ValueError(f"Invalid HiGHS solver option {param!r}={val!r}")

    def relax_model(self):
        """
        Build and solve the continuous relaxation of the current
        model.

        This method creates a HiGHS model in which all integrality
        constraints are relaxed to continuous domains and optimizes
        it. It then clears any cached ``relaxed_model`` property
        so that future accesses re-expose the freshly computed
        object.
        """
        model_copy = self.model.getModel()
        self._relaxed_model = hp.Highs()
        # Options first: HiGHS prints its banner on the first model change if logging is on.
        self._relaxed_model.passOptions(self.model.getOptions())
        self._relaxed_model.passModel(model_copy)

        for i in range(self._relaxed_model.getNumCol()):
            self._relaxed_model.changeColIntegrality(i, hp.HighsVarType.kContinuous)

        self._relaxed_model.optimize()

        self.__dict__.pop("relaxed_model", None)

    def get_best_bound(self) -> float:
        status = self.model.getModelStatus()
        usable = {
            hp.HighsModelStatus.kOptimal,
            hp.HighsModelStatus.kTimeLimit,
            hp.HighsModelStatus.kIterationLimit,
            hp.HighsModelStatus.kSolutionLimit,
            hp.HighsModelStatus.kInterrupt,
            hp.HighsModelStatus.kObjectiveBound,
            hp.HighsModelStatus.kObjectiveTarget,
        }
        # These termination codes are absent from older supported HiGHS releases.
        for name in ("kHighsInterrupt", "kMemoryLimit"):
            value = getattr(hp.HighsModelStatus, name, None)
            if value is not None:
                usable.add(value)
        info = self.model.getInfo()
        if status not in usable or not info.valid:
            raise RuntimeError(f"HiGHS terminated with status {status}; no usable best bound.")
        is_mip = any(t != hp.HighsVarType.kContinuous for t in self.model.getLp().integrality_)
        if is_mip:
            bound = info.mip_dual_bound
        elif status == hp.HighsModelStatus.kOptimal and self.has_solution():
            bound = info.objective_function_value
        else:
            raise RuntimeError(f"HiGHS status {status} has no available LP best bound.")
        if not np.isfinite(bound):
            raise RuntimeError(f"HiGHS status {status} has no finite best bound.")
        return float(bound)

    def has_solution(self) -> bool:
        """Whether the last solve produced a feasible solution."""
        return self.model.getInfo().primal_solution_status == 2

    def relaxed_is_optimal(self) -> bool:
        """Whether the last ``relax_model`` solve ended optimal."""
        return self._relaxed_model.getModelStatus() == hp.HighsModelStatus.kOptimal

    def add_var(self, var: LinearExpression, lb: float = -np.inf, ub: float = np.inf):
        """
        Create and register a list of HiGHS variables corresponding to the
        given ``LinearExpression``.

        This method creates a list of HiGHS variables using ``hp.Highs().addIntegral``
        or ``hp.Highs().addVariable`` to create one or more decision variables whose
        shape matches that of the provided ``LinearExpression``. It also records
        the mapping between the ``LinearExpression`` object and its associated list of
        HiGHS variables for subsequent reference. The created variable(s) have the specified
        bounds and types.

        Parameters
        ----------
        var: LinearExpression
            LinearExpression object to be registered.

        lb: float, default=-np.inf
            Lower bound for all the elements of the matrix variable.

        ub: float, default=np.inf
            Upper bound for all the elements of the matrix variable.

        Returns
        -------
        None
        """

        # HiGHS has no branching priorities, so ``var.branch_priority`` is ignored.
        def _create_single_var(_lb, _ub, is_binary):
            """Create a HiGHS variable using the specified bounds and integrality."""
            self._var_counter += 1
            if is_binary:
                return self.model.addIntegral(lb=0, ub=1, name=f"v{self._var_counter}")
            else:
                return self.model.addVariable(lb=_lb, ub=_ub, name=f"v{self._var_counter}")

        lb = np.full(var.ndim, lb)
        ub = np.full(var.ndim, ub)
        for k, v in var.decomposition_dict.items():
            self.vector_var_map[k] = [
                _create_single_var(lb[i], ub[i], k.is_binary) for i in range(k.ndim)
            ]

    def expression_to_solver(
        self, expression: LinearExpression
    ) -> np.ndarray[hp.highs.highs_linear_expression]:
        r"""
        Convert a ``LinearExpression`` into a list of HiGHS  ``highs_linear_expression``.

        Given a ``LinearExpression`` object with a decomposition of the form
        ``expression.decomposition_dict = {key_1: matrix_1, ..., key_n: matrix_n}``,
        this method constructs the equivalent solver-side representation by multiplying each
        coefficient ``matrix_i`` with its corresponding HiGHS variable and
        summing the results:

        $$
            \sum_{i}^{N} matrix_i @ variable_i
        $$

        The resulting object is a NumPy array of HiGHS linear expressions with the same
        dimensionality as the input ``expression``, the elements of which can be
        used in model constraints or objective functions.

        Parameters
        ----------
        expression: LinearExpression
            The LinearExpression object to be converted into its corresponding HiGHS form.

        Returns
        -------
        np.ndarray[hp.highs.highs_linear_expression]
            An array of HiGHS linear expressions corresponding to the input ``expression``,
            representing its solver-side linear form.
        """
        highs_var = np.zeros(expression.ndim)
        for k, v in expression.decomposition_dict.items():
            v = v.toarray() if spa.issparse(v) else np.array(v)
            if self.vector_var_map.get(k, None) is None:
                self.add_var(k)

            highs_var = highs_var + v @ self.vector_var_map[k]

        return highs_var + expression.const

    def constraints_to_solver(self, constraints: List[Constraint]) -> None:
        """
        Add multiple constraints to the HiGHS model.

        Parameters
        ----------
        constraints: list[Constraint]
            A list of constraint expressions, typically of the form ``lhs == rhs``,
            ``lhs <= rhs``, or ``lhs >= rhs`` where ``lhs`` and ``rhs`` are of type
            ``LinearExpression``.

        Returns
        -------
        None

        """
        for i, constr in enumerate(constraints):
            self.add_constr(constr)

    def update_bounds(
        self, expression: LinearExpression, lb: float | np.ndarray, ub: float | np.ndarray
    ):
        """
        Update the lower and upper bounds of the solver variables corresponding
        to the given ``LinearExpression``.

        Parameters
        ----------
        expression: LinearExpression
            A ``LinearExpression`` object registered with the solver whose
            underlying variable bounds are to be updated.

            lb: float or np.ndarray
            New lower bound(s).

            ub: float or np.ndarray
            New upper bound(s).

        """
        for i, var in enumerate(self.expression_to_solver(expression)):
            var.bounds = (lb[i], ub[i])

    def add_constr(self, constraint: Constraint | hp.highs.highs_linear_expression):
        """
        Add a single constraint to the HiGHS model.

        Parameters
        ----------
        constraint: Constraint | hp.highs.highs_linear_expression
            A compatible HiGHS constraint expression. Typical forms include
            comparisons like ``lhs == rhs``, ``lhs <= rhs``, where ``lhs``,
            and ``rhs`` are of type ``LinearExpression`` or ``hp.highs.highs_linear_expression``.

        Returns
        -------
        None
        """
        try:
            self.model.addConstr(constraint)
        except AttributeError:
            # Expression on the right-hand side of the constraint object is always 0.
            for i, _lhs in enumerate(self.expression_to_solver(constraint.lhs)):
                match constraint.operator:
                    case "=":
                        self.model.addConstr(_lhs == 0)
                    case "<=":
                        self.model.addConstr(_lhs <= 0)

    def set_objective(self, expression: LinearExpression, sense: str = None) -> None:
        """
        Set the optimization objective for the underlying HiGHS model.

        Parameters
        ----------
        expression: LinearExpression or 0
            A ``LinearExpression`` object registered with the solver representing the
            optimization objective, or 0 if setting a zero objective.

        sense: {"minimize", "maximize"}
            Objective sense specified as a string.

        Returns
        -------
        None

        Raises
        ------
        TypeError
            If provided objective ``sense`` is not ``"minimize"`` or ``"maximize"``.

        ValueError
            If provided expression is anything other than a ``LinearExpression`` object or a 0.

        """
        match sense:
            case "minimize":
                sense = ObjSense.kMinimize
            case "maximize":
                sense = ObjSense.kMaximize
            case None:
                sense = None
            case _:
                raise TypeError(f"Objective sense must be either minimize or maximize, got {sense}")

        try:
            if int(expression) != 0:
                raise ValueError(
                    f"Objective function must be either a LinearExpression object or a 0, got {expression}"
                )
            for i in range(self.model.numVariables):
                self.model.changeColCost(i, 0)
        except TypeError:
            self.model.setObjective(self.expression_to_solver(expression)[0], sense=sense)

    # sum_squares (quadratic) constraints are unsupported on HiGHS; the rejecting
    # default in MIPSolver.add_sum_squares_constr applies (it names the solver, so no
    # HiGHS-specific override is needed).

    def optimize(self):
        """
        Optimize the underlying HiGHS model.

        Notes
        -----
        The algorithm used depends on the model objective and
        on the parameter settings
        """
        self.model.optimize()

    def get_solution_value(self, expression: LinearExpression) -> np.ndarray:
        """
        Compute the solution value for the given ``LinearExpression``
        object.

        For each leaf expression in the decomposition of the given
        ``LinearExpression`` object, this method retrieves the corresponding
        solver variables from the model and computes the dot product between
        the coefficient matrix and the vector of variable values associated with
        that leaf.

        Parameters
        ----------
        expression: LinearExpression
            A ``LinearExpression`` object registered with the solver. Each element in
            its decomposition must have a corresponding entry in ``self.vector_var_map``.

        Returns
        -------
        np.ndarray[float]
            A 1-D array of shape ``(expression.ndim,)`` containing the evaluated values
            of the expression at the current model solution.
        """
        if not self.has_solution():
            raise RuntimeError("No feasible incumbent available.")
        out = np.array(expression.const, dtype=float, copy=True)
        for k, v in expression.decomposition_dict.items():
            out += v @ np.array([self.model.val(i) for i in self.vector_var_map[k]])
        return out

    def get_relaxed_variables(
        self, expression: LinearExpression
    ) -> Tuple[list[hp.highs.highs_var], np.ndarray]:
        """
        Reconstruct solver-side variables from a ``LinearExpression`` object
        using the relaxed HiGHS model.

        Parameters
        ----------
        expression: LinearExpression
            A ``LinearExpression`` object registered with the solver. Each element in
            its decomposition must have a corresponding entry in ``self.vector_var_map``.

        Returns
        -------
        Tuple[list[hp.highs.highs_var], np.ndarray]
            A list of HiGHS variables formed by stacking the individual
            ``hp.highs.highs_var`` objects corresponding to each element of
            the input ``LinearExpression``, and a NumPy array containing
            the current value of those stacked variables.
        """
        # If no relaxed model currently exists, build and solve it first.
        if self.relaxed_model is None:
            self.relax_model()

        highs_var, values = [], []
        for k in expression.decomposition_dict.keys():
            _hp_var = self.vector_var_map[k]
            highs_var.extend(_hp_var)
            values.extend(self.relaxed_model.variableValues(_hp_var))

        return highs_var, np.array(values)

    def is_optimal(self) -> bool:
        """
        Check whether the underlying HiGHS model has reached an optimal solution.
        """
        return self.model.getModelStatus() == hp.HighsModelStatus.kOptimal

    def copy(self) -> "HiGHSSolver":
        """
        Create an independent copy of this HiGHS solver.

        Returns a new HiGHSSolver with a copied model and rebuilt
        vector_var_map pointing to variables in the copied model.

        Returns
        -------
        HiGHSSolver
            A new solver instance with copied state.
        """
        new_solver = HiGHSSolver.__new__(HiGHSSolver)
        new_solver._relaxed_model = None
        new_solver.solver_params = self.solver_params
        new_solver.model = hp.Highs()
        # Options first: HiGHS prints its banner on the first model change if logging is on.
        new_solver.model.passOptions(self.model.getOptions())
        new_solver._var_counter = self._var_counter

        # Copy the model data
        model_data = self.model.getModel()
        new_solver.model.passModel(model_data)

        # Build a name-to-variable map for the new model
        new_vars = new_solver.model.getVariables()
        name_to_var = {v.name: v for v in new_vars}

        # Rebuild vector_var_map using variable names
        new_solver.vector_var_map = {}
        for k, var_list in self.vector_var_map.items():
            new_solver.vector_var_map[k] = [name_to_var[v.name] for v in var_list]

        return new_solver

add_constr

add_constr(constraint: Constraint | highs_linear_expression)

Add a single constraint to the HiGHS model.

Parameters:

Name	Type	Description	Default
constraint	Constraint | highs_linear_expression	A compatible HiGHS constraint expression. Typical forms include comparisons like lhs == rhs, lhs <= rhs, where lhs, and rhs are of type LinearExpression or hp.highs.highs_linear_expression.	required

Returns:

Type	Description
None

Source code in fomtools/solvers/highs.py

def add_constr(self, constraint: Constraint | hp.highs.highs_linear_expression):
    """
    Add a single constraint to the HiGHS model.

    Parameters
    ----------
    constraint: Constraint | hp.highs.highs_linear_expression
        A compatible HiGHS constraint expression. Typical forms include
        comparisons like ``lhs == rhs``, ``lhs <= rhs``, where ``lhs``,
        and ``rhs`` are of type ``LinearExpression`` or ``hp.highs.highs_linear_expression``.

    Returns
    -------
    None
    """
    try:
        self.model.addConstr(constraint)
    except AttributeError:
        # Expression on the right-hand side of the constraint object is always 0.
        for i, _lhs in enumerate(self.expression_to_solver(constraint.lhs)):
            match constraint.operator:
                case "=":
                    self.model.addConstr(_lhs == 0)
                case "<=":
                    self.model.addConstr(_lhs <= 0)

add_var

add_var(var: LinearExpression, lb: float = -inf, ub: float = inf)

Create and register a list of HiGHS variables corresponding to the given LinearExpression.

This method creates a list of HiGHS variables using hp.Highs().addIntegral or hp.Highs().addVariable to create one or more decision variables whose shape matches that of the provided LinearExpression. It also records the mapping between the LinearExpression object and its associated list of HiGHS variables for subsequent reference. The created variable(s) have the specified bounds and types.

Parameters:

Name	Type	Description	Default
var	LinearExpression	LinearExpression object to be registered.	required
lb	float	Lower bound for all the elements of the matrix variable.	-inf
ub	float	Upper bound for all the elements of the matrix variable.	inf

Returns:

Type	Description
None

Source code in fomtools/solvers/highs.py

def add_var(self, var: LinearExpression, lb: float = -np.inf, ub: float = np.inf):
    """
    Create and register a list of HiGHS variables corresponding to the
    given ``LinearExpression``.

    This method creates a list of HiGHS variables using ``hp.Highs().addIntegral``
    or ``hp.Highs().addVariable`` to create one or more decision variables whose
    shape matches that of the provided ``LinearExpression``. It also records
    the mapping between the ``LinearExpression`` object and its associated list of
    HiGHS variables for subsequent reference. The created variable(s) have the specified
    bounds and types.

    Parameters
    ----------
    var: LinearExpression
        LinearExpression object to be registered.

    lb: float, default=-np.inf
        Lower bound for all the elements of the matrix variable.

    ub: float, default=np.inf
        Upper bound for all the elements of the matrix variable.

    Returns
    -------
    None
    """

    # HiGHS has no branching priorities, so ``var.branch_priority`` is ignored.
    def _create_single_var(_lb, _ub, is_binary):
        """Create a HiGHS variable using the specified bounds and integrality."""
        self._var_counter += 1
        if is_binary:
            return self.model.addIntegral(lb=0, ub=1, name=f"v{self._var_counter}")
        else:
            return self.model.addVariable(lb=_lb, ub=_ub, name=f"v{self._var_counter}")

    lb = np.full(var.ndim, lb)
    ub = np.full(var.ndim, ub)
    for k, v in var.decomposition_dict.items():
        self.vector_var_map[k] = [
            _create_single_var(lb[i], ub[i], k.is_binary) for i in range(k.ndim)
        ]

constraints_to_solver

constraints_to_solver(constraints: List[Constraint]) -> None

Add multiple constraints to the HiGHS model.

Parameters:

Name	Type	Description	Default
constraints	List[Constraint]	A list of constraint expressions, typically of the form lhs == rhs, lhs <= rhs, or lhs >= rhs where lhs and rhs are of type LinearExpression.	required

Returns:

Type	Description
None

Source code in fomtools/solvers/highs.py

def constraints_to_solver(self, constraints: List[Constraint]) -> None:
    """
    Add multiple constraints to the HiGHS model.

    Parameters
    ----------
    constraints: list[Constraint]
        A list of constraint expressions, typically of the form ``lhs == rhs``,
        ``lhs <= rhs``, or ``lhs >= rhs`` where ``lhs`` and ``rhs`` are of type
        ``LinearExpression``.

    Returns
    -------
    None

    """
    for i, constr in enumerate(constraints):
        self.add_constr(constr)

copy

copy() -> HiGHSSolver

Create an independent copy of this HiGHS solver.

Returns a new HiGHSSolver with a copied model and rebuilt vector_var_map pointing to variables in the copied model.

Returns:

Type	Description
HiGHSSolver	A new solver instance with copied state.

Source code in fomtools/solvers/highs.py

def copy(self) -> "HiGHSSolver":
    """
    Create an independent copy of this HiGHS solver.

    Returns a new HiGHSSolver with a copied model and rebuilt
    vector_var_map pointing to variables in the copied model.

    Returns
    -------
    HiGHSSolver
        A new solver instance with copied state.
    """
    new_solver = HiGHSSolver.__new__(HiGHSSolver)
    new_solver._relaxed_model = None
    new_solver.solver_params = self.solver_params
    new_solver.model = hp.Highs()
    # Options first: HiGHS prints its banner on the first model change if logging is on.
    new_solver.model.passOptions(self.model.getOptions())
    new_solver._var_counter = self._var_counter

    # Copy the model data
    model_data = self.model.getModel()
    new_solver.model.passModel(model_data)

    # Build a name-to-variable map for the new model
    new_vars = new_solver.model.getVariables()
    name_to_var = {v.name: v for v in new_vars}

    # Rebuild vector_var_map using variable names
    new_solver.vector_var_map = {}
    for k, var_list in self.vector_var_map.items():
        new_solver.vector_var_map[k] = [name_to_var[v.name] for v in var_list]

    return new_solver

expression_to_solver

expression_to_solver(expression: LinearExpression) -> ndarray[highs_linear_expression]

Convert a LinearExpression into a list of HiGHS highs_linear_expression.

Given a LinearExpression object with a decomposition of the form expression.decomposition_dict = {key_1: matrix_1, ..., key_n: matrix_n}, this method constructs the equivalent solver-side representation by multiplying each coefficient matrix_i with its corresponding HiGHS variable and summing the results:

\[ \sum_{i}^{N} matrix_i @ variable_i \]

The resulting object is a NumPy array of HiGHS linear expressions with the same dimensionality as the input expression, the elements of which can be used in model constraints or objective functions.

Parameters:

Name	Type	Description	Default
expression	LinearExpression	The LinearExpression object to be converted into its corresponding HiGHS form.	required

Returns:

Type	Description
ndarray[highs_linear_expression]	An array of HiGHS linear expressions corresponding to the input expression, representing its solver-side linear form.

Source code in fomtools/solvers/highs.py

def expression_to_solver(
    self, expression: LinearExpression
) -> np.ndarray[hp.highs.highs_linear_expression]:
    r"""
    Convert a ``LinearExpression`` into a list of HiGHS  ``highs_linear_expression``.

    Given a ``LinearExpression`` object with a decomposition of the form
    ``expression.decomposition_dict = {key_1: matrix_1, ..., key_n: matrix_n}``,
    this method constructs the equivalent solver-side representation by multiplying each
    coefficient ``matrix_i`` with its corresponding HiGHS variable and
    summing the results:

    $$
        \sum_{i}^{N} matrix_i @ variable_i
    $$

    The resulting object is a NumPy array of HiGHS linear expressions with the same
    dimensionality as the input ``expression``, the elements of which can be
    used in model constraints or objective functions.

    Parameters
    ----------
    expression: LinearExpression
        The LinearExpression object to be converted into its corresponding HiGHS form.

    Returns
    -------
    np.ndarray[hp.highs.highs_linear_expression]
        An array of HiGHS linear expressions corresponding to the input ``expression``,
        representing its solver-side linear form.
    """
    highs_var = np.zeros(expression.ndim)
    for k, v in expression.decomposition_dict.items():
        v = v.toarray() if spa.issparse(v) else np.array(v)
        if self.vector_var_map.get(k, None) is None:
            self.add_var(k)

        highs_var = highs_var + v @ self.vector_var_map[k]

    return highs_var + expression.const

get_relaxed_variables

get_relaxed_variables(expression: LinearExpression) -> Tuple[list[highs_var], ndarray]

Reconstruct solver-side variables from a LinearExpression object using the relaxed HiGHS model.

Parameters:

Name	Type	Description	Default
expression	LinearExpression	A LinearExpression object registered with the solver. Each element in its decomposition must have a corresponding entry in self.vector_var_map.	required

Returns:

Type	Description
Tuple[list[highs_var], ndarray]	A list of HiGHS variables formed by stacking the individual hp.highs.highs_var objects corresponding to each element of the input LinearExpression, and a NumPy array containing the current value of those stacked variables.

Source code in fomtools/solvers/highs.py

def get_relaxed_variables(
    self, expression: LinearExpression
) -> Tuple[list[hp.highs.highs_var], np.ndarray]:
    """
    Reconstruct solver-side variables from a ``LinearExpression`` object
    using the relaxed HiGHS model.

    Parameters
    ----------
    expression: LinearExpression
        A ``LinearExpression`` object registered with the solver. Each element in
        its decomposition must have a corresponding entry in ``self.vector_var_map``.

    Returns
    -------
    Tuple[list[hp.highs.highs_var], np.ndarray]
        A list of HiGHS variables formed by stacking the individual
        ``hp.highs.highs_var`` objects corresponding to each element of
        the input ``LinearExpression``, and a NumPy array containing
        the current value of those stacked variables.
    """
    # If no relaxed model currently exists, build and solve it first.
    if self.relaxed_model is None:
        self.relax_model()

    highs_var, values = [], []
    for k in expression.decomposition_dict.keys():
        _hp_var = self.vector_var_map[k]
        highs_var.extend(_hp_var)
        values.extend(self.relaxed_model.variableValues(_hp_var))

    return highs_var, np.array(values)

get_solution_value

get_solution_value(expression: LinearExpression) -> ndarray

Compute the solution value for the given LinearExpression object.

For each leaf expression in the decomposition of the given LinearExpression object, this method retrieves the corresponding solver variables from the model and computes the dot product between the coefficient matrix and the vector of variable values associated with that leaf.

Parameters:

Name	Type	Description	Default
expression	LinearExpression	A LinearExpression object registered with the solver. Each element in its decomposition must have a corresponding entry in self.vector_var_map.	required

Returns:

Type	Description
ndarray[float]	A 1-D array of shape (expression.ndim,) containing the evaluated values of the expression at the current model solution.

Source code in fomtools/solvers/highs.py

def get_solution_value(self, expression: LinearExpression) -> np.ndarray:
    """
    Compute the solution value for the given ``LinearExpression``
    object.

    For each leaf expression in the decomposition of the given
    ``LinearExpression`` object, this method retrieves the corresponding
    solver variables from the model and computes the dot product between
    the coefficient matrix and the vector of variable values associated with
    that leaf.

    Parameters
    ----------
    expression: LinearExpression
        A ``LinearExpression`` object registered with the solver. Each element in
        its decomposition must have a corresponding entry in ``self.vector_var_map``.

    Returns
    -------
    np.ndarray[float]
        A 1-D array of shape ``(expression.ndim,)`` containing the evaluated values
        of the expression at the current model solution.
    """
    if not self.has_solution():
        raise RuntimeError("No feasible incumbent available.")
    out = np.array(expression.const, dtype=float, copy=True)
    for k, v in expression.decomposition_dict.items():
        out += v @ np.array([self.model.val(i) for i in self.vector_var_map[k]])
    return out

has_solution

has_solution() -> bool

Whether the last solve produced a feasible solution.

Source code in fomtools/solvers/highs.py

def has_solution(self) -> bool:
    """Whether the last solve produced a feasible solution."""
    return self.model.getInfo().primal_solution_status == 2

is_optimal

is_optimal() -> bool

Check whether the underlying HiGHS model has reached an optimal solution.

Source code in fomtools/solvers/highs.py

def is_optimal(self) -> bool:
    """
    Check whether the underlying HiGHS model has reached an optimal solution.
    """
    return self.model.getModelStatus() == hp.HighsModelStatus.kOptimal

optimize

optimize()

Optimize the underlying HiGHS model.

Notes

The algorithm used depends on the model objective and on the parameter settings

Source code in fomtools/solvers/highs.py

def optimize(self):
    """
    Optimize the underlying HiGHS model.

    Notes
    -----
    The algorithm used depends on the model objective and
    on the parameter settings
    """
    self.model.optimize()

relax_model

relax_model()

Build and solve the continuous relaxation of the current model.

This method creates a HiGHS model in which all integrality constraints are relaxed to continuous domains and optimizes it. It then clears any cached relaxed_model property so that future accesses re-expose the freshly computed object.

Source code in fomtools/solvers/highs.py

def relax_model(self):
    """
    Build and solve the continuous relaxation of the current
    model.

    This method creates a HiGHS model in which all integrality
    constraints are relaxed to continuous domains and optimizes
    it. It then clears any cached ``relaxed_model`` property
    so that future accesses re-expose the freshly computed
    object.
    """
    model_copy = self.model.getModel()
    self._relaxed_model = hp.Highs()
    # Options first: HiGHS prints its banner on the first model change if logging is on.
    self._relaxed_model.passOptions(self.model.getOptions())
    self._relaxed_model.passModel(model_copy)

    for i in range(self._relaxed_model.getNumCol()):
        self._relaxed_model.changeColIntegrality(i, hp.HighsVarType.kContinuous)

    self._relaxed_model.optimize()

    self.__dict__.pop("relaxed_model", None)

relaxed_is_optimal

relaxed_is_optimal() -> bool

Whether the last relax_model solve ended optimal.

Source code in fomtools/solvers/highs.py

def relaxed_is_optimal(self) -> bool:
    """Whether the last ``relax_model`` solve ended optimal."""
    return self._relaxed_model.getModelStatus() == hp.HighsModelStatus.kOptimal

set_objective

set_objective(expression: LinearExpression, sense: str = None) -> None

Set the optimization objective for the underlying HiGHS model.

Parameters:

Name	Type	Description	Default
expression	LinearExpression	A LinearExpression object registered with the solver representing the optimization objective, or 0 if setting a zero objective.	required
sense	str	Objective sense specified as a string.	None

Returns:

Type	Description
None

Raises:

Type	Description
TypeError	If provided objective sense is not "minimize" or "maximize".
ValueError	If provided expression is anything other than a LinearExpression object or a 0.

Source code in fomtools/solvers/highs.py

def set_objective(self, expression: LinearExpression, sense: str = None) -> None:
    """
    Set the optimization objective for the underlying HiGHS model.

    Parameters
    ----------
    expression: LinearExpression or 0
        A ``LinearExpression`` object registered with the solver representing the
        optimization objective, or 0 if setting a zero objective.

    sense: {"minimize", "maximize"}
        Objective sense specified as a string.

    Returns
    -------
    None

    Raises
    ------
    TypeError
        If provided objective ``sense`` is not ``"minimize"`` or ``"maximize"``.

    ValueError
        If provided expression is anything other than a ``LinearExpression`` object or a 0.

    """
    match sense:
        case "minimize":
            sense = ObjSense.kMinimize
        case "maximize":
            sense = ObjSense.kMaximize
        case None:
            sense = None
        case _:
            raise TypeError(f"Objective sense must be either minimize or maximize, got {sense}")

    try:
        if int(expression) != 0:
            raise ValueError(
                f"Objective function must be either a LinearExpression object or a 0, got {expression}"
            )
        for i in range(self.model.numVariables):
            self.model.changeColCost(i, 0)
    except TypeError:
        self.model.setObjective(self.expression_to_solver(expression)[0], sense=sense)

set_param

set_param(param: str, val)

Configure a solver option, passed verbatim to HiGHS.

The option name must be a native HiGHS option (e.g. time_limit, mip_rel_gap, output_flag); no name translation is performed.

Parameters:

Name	Type	Description	Default
param	str	The native HiGHS option name to set.	required
val		The value to assign to the specified option.	required

Raises:

Type	Description
ValueError	If HiGHS rejects the option name or value.

Source code in fomtools/solvers/highs.py

def set_param(self, param: str, val):
    """
    Configure a solver option, passed verbatim to HiGHS.

    The option name must be a native HiGHS option (e.g. ``time_limit``,
    ``mip_rel_gap``, ``output_flag``); no name translation is performed.

    Parameters
    ----------
    param: str
        The native HiGHS option name to set.
    val:
        The value to assign to the specified option.

    Raises
    ------
    ValueError
        If HiGHS rejects the option name or value.
    """
    status = self.model.setOptionValue(param, val)
    if status == hp.HighsStatus.kError:
        raise ValueError(f"Invalid HiGHS solver option {param!r}={val!r}")

update_bounds

update_bounds(expression: LinearExpression, lb: float | ndarray, ub: float | ndarray)

Update the lower and upper bounds of the solver variables corresponding to the given LinearExpression.

Parameters:

Name	Type	Description	Default
expression	LinearExpression	A LinearExpression object registered with the solver whose underlying variable bounds are to be updated. lb: float or np.ndarray New lower bound(s). ub: float or np.ndarray New upper bound(s).	required

Source code in fomtools/solvers/highs.py

def update_bounds(
    self, expression: LinearExpression, lb: float | np.ndarray, ub: float | np.ndarray
):
    """
    Update the lower and upper bounds of the solver variables corresponding
    to the given ``LinearExpression``.

    Parameters
    ----------
    expression: LinearExpression
        A ``LinearExpression`` object registered with the solver whose
        underlying variable bounds are to be updated.

        lb: float or np.ndarray
        New lower bound(s).

        ub: float or np.ndarray
        New upper bound(s).

    """
    for i, var in enumerate(self.expression_to_solver(expression)):
        var.bounds = (lb[i], ub[i])

ClarabelSolver

Bases: SDPSolver

Solver wrapper for the Clarabel conic optimization solver.

This class provides an interface to solve semidefinite programs (SDPs) using the Clarabel solver. It converts the SDP problem data (equality constraints, inequality constraints, and PSD cone constraints) into Clarabel's native format and invokes the solver.

The SDP is formulated as: maximize tr(C @ X) subject to tr(Aeq[i] @ X) = beq[i] for all i tr(Aineq[j] @ X) <= bineq[j] for all j X[cone_indices] >> 0 for each PSD cone

Attributes:

Name	Type	Description
q	ndarray	Vectorized objective matrix.
Aeq	csc_matrix	Stacked equality constraint matrices (vectorized).
beq	ndarray	Equality constraint right-hand sides.
Aineq	csc_matrix	Stacked inequality constraint matrices (vectorized).
bineq	ndarray	Inequality constraint right-hand sides.
M	csc_matrix	PSD cone selection matrices.
cone_dims	dict	Cone dimensions with keys 'zero', 'nonneg', 'psd'.

Source code in fomtools/solvers/clarabel.py

class ClarabelSolver(SDPSolver):
    """
    Solver wrapper for the Clarabel conic optimization solver.

    This class provides an interface to solve semidefinite programs (SDPs) using
    the Clarabel solver. It converts the SDP problem data (equality constraints,
    inequality constraints, and PSD cone constraints) into Clarabel's native format
    and invokes the solver.

    The SDP is formulated as:
        maximize    tr(C @ X)
        subject to  tr(Aeq[i] @ X) = beq[i]      for all i
                    tr(Aineq[j] @ X) <= bineq[j]  for all j
                    X[cone_indices] >> 0          for each PSD cone

    Attributes
    ----------
    q : np.ndarray
        Vectorized objective matrix.
    Aeq : scipy.sparse.csc_matrix
        Stacked equality constraint matrices (vectorized).
    beq : np.ndarray
        Equality constraint right-hand sides.
    Aineq : scipy.sparse.csc_matrix
        Stacked inequality constraint matrices (vectorized).
    bineq : np.ndarray
        Inequality constraint right-hand sides.
    M : scipy.sparse.csc_matrix
        PSD cone selection matrices.
    cone_dims : dict
        Cone dimensions with keys 'zero', 'nonneg', 'psd'.
    """

    def __init__(self, solver_params=None):
        """
        Initialize the Clarabel solver.

        Parameters
        ----------
        solver_params : dict, optional
            Dictionary of native Clarabel settings to configure. Stored
            and applied to the ``DefaultSettings`` object at solve time.
        """
        self.solver_params = solver_params or {}

    def setup_clarabel_data(
        self,
        mat_problem_dim,
        C,
        Aeq_mats,
        beq_consts,
        Aineq_mats,
        bineq_consts,
        psd_cone_indices,
    ):
        """
        Convert SDP problem data to Clarabel's native format.

        This method transforms the matrix-form SDP specification into the
        vectorized conic form expected by Clarabel. Symmetric matrices are
        vectorized using column-major upper triangular format with appropriate
        scaling for off-diagonal elements.

        Parameters
        ----------
        mat_problem_dim : int
            Dimension of the symmetric matrix variable X.
        C : scipy.sparse.csc_matrix
            Objective matrix for maximizing tr(C @ X).
        Aeq_mats : list[scipy.sparse.csc_matrix]
            List of equality constraint matrices A_i for tr(A_i @ X) = b_i.
        beq_consts : list[float]
            Right-hand sides for equality constraints.
        Aineq_mats : list[scipy.sparse.csc_matrix]
            List of inequality constraint matrices for tr(A_j @ X) <= b_j.
        bineq_consts : list[float]
            Right-hand sides for inequality constraints.
        psd_cone_indices : list[set]
            List of index sets defining PSD cone constraints on submatrices.
        """
        vec_problem_dim = mat_problem_dim * (mat_problem_dim + 1) // 2

        q = -sparse_symmetric_vectorize(C, scale_factor=2)
        # TODO: see if this can be replaced by a sparse 1 by n vector instead
        q = (
            q.todense()
            .getA()
            .reshape(
                -1,
            )
        )
        if q.shape[0] != vec_problem_dim:
            raise ValueError("Objective dimension does not match the SDP variable.")

        symmetrized_Aeq_mats = [sparse_symmetric_vectorize(A, scale_factor=2) for A in Aeq_mats]
        Aeq = spa.vstack(symmetrized_Aeq_mats)
        beq = np.array(beq_consts)

        if Aeq.shape != (beq.shape[0], vec_problem_dim):
            raise ValueError("Equality data have inconsistent dimensions.")

        symmetrized_Aineq_mats = [sparse_symmetric_vectorize(A, scale_factor=2) for A in Aineq_mats]
        Aineq = spa.vstack(symmetrized_Aineq_mats)
        bineq = np.array(bineq_consts)

        if Aineq.shape != (bineq.shape[0], vec_problem_dim):
            raise ValueError("Inequality data have inconsistent dimensions.")

        psd_cone_sizes = []
        psd_cone_selectors = []

        for cone in psd_cone_indices:
            n = len(cone)
            psd_cone_sizes.append(n)
            M = -get_submatrix_selection_matrix(mat_problem_dim, list(cone))
            psd_cone_selectors.append(M)

        M = spa.vstack(psd_cone_selectors)

        cone_dims = {
            "zero": Aeq.shape[0],
            "nonneg": Aineq.shape[0],
            "psd": psd_cone_sizes,
        }

        self.q = q
        self.Aeq = Aeq
        self.beq = beq
        self.Aineq = Aineq
        self.bineq = bineq
        self.M = M
        self.cone_dims = cone_dims

    def solve(self, **solve_kwargs):
        """
        Solve the SDP and return its accepted dual upper bound.

        Assembles the conic problem from the prepared data and invokes
        the Clarabel solver. The problem is formulated as a minimization
        in Clarabel's standard form, so the returned value is negated
        to give a dual upper bound for the maximization problem.
        ``Solved`` is accepted; ``AlmostSolved`` is accepted with a warning, since
        the dual objective then holds only to Clarabel's reduced tolerances. Any
        other termination (limits, infeasibility, numerical errors) raises.

        Parameters
        ----------
        **solve_kwargs
            Native Clarabel settings applied verbatim to the
            ``DefaultSettings`` object (e.g. ``max_iter``,
            ``tol_gap_abs``). These override matching entries in
            ``solver_params`` for this solve.

        Returns
        -------
        float
            The finite dual objective at the requested numerical tolerances.

        Raises
        ------
        ValueError
            If Clarabel rejects a settings name or value.
        """
        q = self.q
        b_psd = [n * (n + 1) // 2 for n in self.cone_dims["psd"]]
        b = np.concatenate([self.beq, self.bineq, np.zeros(sum(b_psd))])
        A = spa.vstack([self.Aeq, self.Aineq, self.M])
        if A.shape != (b.shape[0], q.shape[0]):
            raise ValueError("Conic data have inconsistent dimensions.")

        cones = [
            clarabel.ZeroConeT(self.beq.shape[0]),
            clarabel.NonnegativeConeT(self.bineq.shape[0]),
        ]
        cones += [clarabel.PSDTriangleConeT(n) for n in self.cone_dims["psd"]]

        xdim = q.shape[0]
        P = spa.csc_matrix((xdim, xdim))
        settings = clarabel.DefaultSettings()
        settings.verbose = False

        # Apply native Clarabel settings; solve-time kwargs override construction params.
        merged = {**self.solver_params, **solve_kwargs}
        for k, v in merged.items():
            try:
                setattr(settings, k, v)
            except (AttributeError, TypeError) as e:
                raise ValueError(f"Invalid Clarabel solver option {k!r}={v!r}: {e}") from e

        solver = clarabel.DefaultSolver(P, q, A, b, cones, settings)
        self.solution = sol = solver.solve()
        if sol.status == clarabel.SolverStatus.AlmostSolved:
            warnings.warn(
                "Clarabel terminated AlmostSolved; the SDP bound holds only to its "
                "reduced tolerances.",
                stacklevel=2,
            )
        elif sol.status != clarabel.SolverStatus.Solved:
            raise RuntimeError(f"Clarabel terminated with {sol.status}; no accepted SDP bound.")
        bound = -float(sol.obj_val_dual)
        if not np.isfinite(bound):
            raise RuntimeError("Clarabel returned a nonfinite SDP dual objective.")
        return bound

__init__

__init__(solver_params=None)

Parameters:

Name	Type	Description	Default
solver_params	dict	Dictionary of native Clarabel settings to configure. Stored and applied to the DefaultSettings object at solve time.	None

Source code in fomtools/solvers/clarabel.py

def __init__(self, solver_params=None):
    """
    Initialize the Clarabel solver.

    Parameters
    ----------
    solver_params : dict, optional
        Dictionary of native Clarabel settings to configure. Stored
        and applied to the ``DefaultSettings`` object at solve time.
    """
    self.solver_params = solver_params or {}

setup_clarabel_data

setup_clarabel_data(mat_problem_dim, C, Aeq_mats, beq_consts, Aineq_mats, bineq_consts, psd_cone_indices)

Convert SDP problem data to Clarabel's native format.

This method transforms the matrix-form SDP specification into the vectorized conic form expected by Clarabel. Symmetric matrices are vectorized using column-major upper triangular format with appropriate scaling for off-diagonal elements.

Parameters:

Name	Type	Description	Default
mat_problem_dim	int	Dimension of the symmetric matrix variable X.	required
C	csc_matrix	Objective matrix for maximizing tr(C @ X).	required
Aeq_mats	list[csc_matrix]	List of equality constraint matrices A_i for tr(A_i @ X) = b_i.	required
beq_consts	list[float]	Right-hand sides for equality constraints.	required
Aineq_mats	list[csc_matrix]	List of inequality constraint matrices for tr(A_j @ X) <= b_j.	required
bineq_consts	list[float]	Right-hand sides for inequality constraints.	required
psd_cone_indices	list[set]	List of index sets defining PSD cone constraints on submatrices.	required

Source code in fomtools/solvers/clarabel.py

def setup_clarabel_data(
    self,
    mat_problem_dim,
    C,
    Aeq_mats,
    beq_consts,
    Aineq_mats,
    bineq_consts,
    psd_cone_indices,
):
    """
    Convert SDP problem data to Clarabel's native format.

    This method transforms the matrix-form SDP specification into the
    vectorized conic form expected by Clarabel. Symmetric matrices are
    vectorized using column-major upper triangular format with appropriate
    scaling for off-diagonal elements.

    Parameters
    ----------
    mat_problem_dim : int
        Dimension of the symmetric matrix variable X.
    C : scipy.sparse.csc_matrix
        Objective matrix for maximizing tr(C @ X).
    Aeq_mats : list[scipy.sparse.csc_matrix]
        List of equality constraint matrices A_i for tr(A_i @ X) = b_i.
    beq_consts : list[float]
        Right-hand sides for equality constraints.
    Aineq_mats : list[scipy.sparse.csc_matrix]
        List of inequality constraint matrices for tr(A_j @ X) <= b_j.
    bineq_consts : list[float]
        Right-hand sides for inequality constraints.
    psd_cone_indices : list[set]
        List of index sets defining PSD cone constraints on submatrices.
    """
    vec_problem_dim = mat_problem_dim * (mat_problem_dim + 1) // 2

    q = -sparse_symmetric_vectorize(C, scale_factor=2)
    # TODO: see if this can be replaced by a sparse 1 by n vector instead
    q = (
        q.todense()
        .getA()
        .reshape(
            -1,
        )
    )
    if q.shape[0] != vec_problem_dim:
        raise ValueError("Objective dimension does not match the SDP variable.")

    symmetrized_Aeq_mats = [sparse_symmetric_vectorize(A, scale_factor=2) for A in Aeq_mats]
    Aeq = spa.vstack(symmetrized_Aeq_mats)
    beq = np.array(beq_consts)

    if Aeq.shape != (beq.shape[0], vec_problem_dim):
        raise ValueError("Equality data have inconsistent dimensions.")

    symmetrized_Aineq_mats = [sparse_symmetric_vectorize(A, scale_factor=2) for A in Aineq_mats]
    Aineq = spa.vstack(symmetrized_Aineq_mats)
    bineq = np.array(bineq_consts)

    if Aineq.shape != (bineq.shape[0], vec_problem_dim):
        raise ValueError("Inequality data have inconsistent dimensions.")

    psd_cone_sizes = []
    psd_cone_selectors = []

    for cone in psd_cone_indices:
        n = len(cone)
        psd_cone_sizes.append(n)
        M = -get_submatrix_selection_matrix(mat_problem_dim, list(cone))
        psd_cone_selectors.append(M)

    M = spa.vstack(psd_cone_selectors)

    cone_dims = {
        "zero": Aeq.shape[0],
        "nonneg": Aineq.shape[0],
        "psd": psd_cone_sizes,
    }

    self.q = q
    self.Aeq = Aeq
    self.beq = beq
    self.Aineq = Aineq
    self.bineq = bineq
    self.M = M
    self.cone_dims = cone_dims

solve

solve(**solve_kwargs)

Solve the SDP and return its accepted dual upper bound.

Assembles the conic problem from the prepared data and invokes the Clarabel solver. The problem is formulated as a minimization in Clarabel's standard form, so the returned value is negated to give a dual upper bound for the maximization problem. Solved is accepted; AlmostSolved is accepted with a warning, since the dual objective then holds only to Clarabel's reduced tolerances. Any other termination (limits, infeasibility, numerical errors) raises.

Parameters:

Name	Type	Description	Default
**solve_kwargs		Native Clarabel settings applied verbatim to the DefaultSettings object (e.g. max_iter, tol_gap_abs). These override matching entries in solver_params for this solve.	{}

Returns:

Type	Description
float	The finite dual objective at the requested numerical tolerances.

Raises:

Type	Description
ValueError	If Clarabel rejects a settings name or value.

Source code in fomtools/solvers/clarabel.py

def solve(self, **solve_kwargs):
    """
    Solve the SDP and return its accepted dual upper bound.

    Assembles the conic problem from the prepared data and invokes
    the Clarabel solver. The problem is formulated as a minimization
    in Clarabel's standard form, so the returned value is negated
    to give a dual upper bound for the maximization problem.
    ``Solved`` is accepted; ``AlmostSolved`` is accepted with a warning, since
    the dual objective then holds only to Clarabel's reduced tolerances. Any
    other termination (limits, infeasibility, numerical errors) raises.

    Parameters
    ----------
    **solve_kwargs
        Native Clarabel settings applied verbatim to the
        ``DefaultSettings`` object (e.g. ``max_iter``,
        ``tol_gap_abs``). These override matching entries in
        ``solver_params`` for this solve.

    Returns
    -------
    float
        The finite dual objective at the requested numerical tolerances.

    Raises
    ------
    ValueError
        If Clarabel rejects a settings name or value.
    """
    q = self.q
    b_psd = [n * (n + 1) // 2 for n in self.cone_dims["psd"]]
    b = np.concatenate([self.beq, self.bineq, np.zeros(sum(b_psd))])
    A = spa.vstack([self.Aeq, self.Aineq, self.M])
    if A.shape != (b.shape[0], q.shape[0]):
        raise ValueError("Conic data have inconsistent dimensions.")

    cones = [
        clarabel.ZeroConeT(self.beq.shape[0]),
        clarabel.NonnegativeConeT(self.bineq.shape[0]),
    ]
    cones += [clarabel.PSDTriangleConeT(n) for n in self.cone_dims["psd"]]

    xdim = q.shape[0]
    P = spa.csc_matrix((xdim, xdim))
    settings = clarabel.DefaultSettings()
    settings.verbose = False

    # Apply native Clarabel settings; solve-time kwargs override construction params.
    merged = {**self.solver_params, **solve_kwargs}
    for k, v in merged.items():
        try:
            setattr(settings, k, v)
        except (AttributeError, TypeError) as e:
            raise ValueError(f"Invalid Clarabel solver option {k!r}={v!r}: {e}") from e

    solver = clarabel.DefaultSolver(P, q, A, b, cones, settings)
    self.solution = sol = solver.solve()
    if sol.status == clarabel.SolverStatus.AlmostSolved:
        warnings.warn(
            "Clarabel terminated AlmostSolved; the SDP bound holds only to its "
            "reduced tolerances.",
            stacklevel=2,
        )
    elif sol.status != clarabel.SolverStatus.Solved:
        raise RuntimeError(f"Clarabel terminated with {sol.status}; no accepted SDP bound.")
    bound = -float(sol.obj_val_dual)
    if not np.isfinite(bound):
        raise RuntimeError("Clarabel returned a nonfinite SDP dual objective.")
    return bound

CuClarabelSolver

Bases: ClarabelSolver

GPU-accelerated solver wrapper for Clarabel using CUDA.

This class extends ClarabelSolver to leverage GPU acceleration via CuPy and Julia's CUDA-enabled Clarabel implementation. It transfers problem data to GPU memory and solves the SDP using cuDSS for sparse linear algebra.

Requirements

• CuPy with CUDA support
• Julia with Clarabel.jl, CUDA.jl packages
• juliacall Python package

Notes

Problem setup uses the parent class's setup_clarabel_data method. Only the solve method differs, utilizing GPU kernels.

Source code in fomtools/solvers/cuclarabel.py

class CuClarabelSolver(ClarabelSolver):
    """
    GPU-accelerated solver wrapper for Clarabel using CUDA.

    This class extends ClarabelSolver to leverage GPU acceleration via CuPy and
    Julia's CUDA-enabled Clarabel implementation. It transfers problem data to
    GPU memory and solves the SDP using cuDSS for sparse linear algebra.

    Requirements
    ------------
    - CuPy with CUDA support
    - Julia with Clarabel.jl, CUDA.jl packages
    - juliacall Python package

    Notes
    -----
    Problem setup uses the parent class's ``setup_clarabel_data`` method.
    Only the ``solve`` method differs, utilizing GPU kernels.
    """

    def __init__(self, solver_params=None):
        """
        Initialize the GPU-accelerated Clarabel solver.

        Parameters
        ----------
        solver_params : dict, optional
            Dictionary of solver-specific parameters to configure.

        Raises
        ------
        ValueError
            If CUDA dependencies (cupy) are not available.
        """
        _check_cuda_available()
        super(CuClarabelSolver, self).__init__(solver_params=solver_params)

    def solve(self, **solve_kwargs):
        """
        Solve the SDP on GPU and return the objective value.

        Transfers the problem data to GPU memory using CuPy and invokes
        Julia's CUDA-enabled Clarabel solver via juliacall. Uses cuDSS
        for the direct linear system solves.

        Parameters
        ----------
        **solve_kwargs
            Accepted for API parity, but native Clarabel settings are
            not yet wired into the Julia GPU path. Passing any option
            raises ``ValueError``.

        Returns
        -------
        float
            The optimal objective value of the SDP maximization problem.

        Notes
        -----
        Requires CUDA-capable GPU and properly configured Julia environment
        with CUDA.jl and Clarabel.jl packages.
        """
        merged = {**self.solver_params, **solve_kwargs}
        if merged:
            raise ValueError(
                "cuClarabel does not support custom solver options yet; "
                f"got {merged}. (Construction-time params are rejected up front "
                "by the SDP canonicalizer.)"
            )
        q = self.q
        b_psd = [n * (n + 1) // 2 for n in self.cone_dims["psd"]]
        b = np.concatenate([self.beq, self.bineq, np.zeros(sum(b_psd))])
        A = spa.vstack([self.Aeq, self.Aineq, self.M])

        ndim = q.shape[0]

        import cupy
        from cupyx.scipy.sparse import csr_matrix as cucsr_matrix
        from juliacall import Main as jl

        jl.seval("using CUDA")
        jl.seval("using Clarabel")
        jl.seval("using LinearAlgebra, SparseArrays, CUDA.CUSPARSE")

        # 2. Move Data to GPU
        qgpu = cupy.array(q)
        Agpu = cucsr_matrix(A)
        bgpu = cupy.array(b)

        jl.seval("global q; global A; global b")

        # 5. Transfer Pointers to Julia
        # Note: We initialize P as an empty GPU matrix since it is Zero
        jl.seval(f"""
            P = CuSparseMatrixCSR(sparse(Float64[], Float64[], Float64[], {ndim}, {ndim}))
        """)

        jl.q = jl.Clarabel.cupy_to_cuvector(jl.Float64, int(qgpu.data.ptr), qgpu.size)

        jl.A = jl.Clarabel.cupy_to_cucsrmat(
            jl.Float64,
            int(Agpu.data.data.ptr),
            int(Agpu.indices.data.ptr),
            int(Agpu.indptr.data.ptr),
            *Agpu.shape,
            Agpu.nnz,
        )
        jl.b = jl.Clarabel.cupy_to_cuvector(jl.Float64, int(bgpu.data.ptr), bgpu.size)

        # 6. Setup Cones
        jl.seval("cones = Clarabel.SupportedCone[]")
        jl.push_b(jl.cones, jl.Clarabel.ZeroConeT(self.beq.shape[0]))
        jl.push_b(jl.cones, jl.Clarabel.NonnegativeConeT(self.bineq.shape[0]))
        for n in self.cone_dims["psd"]:
            jl.push_b(jl.cones, jl.Clarabel.PSDTriangleConeT(n))

        # 7. Solve
        # The 'MethodError' will vanish because the Extension now provides the GPU signature.
        # settings = Clarabel.Settings(direct_solve_method = :cudss, verbose = true)
        jl.seval("""
            settings = Clarabel.Settings(direct_solve_method = :cudss, verbose = false)
            solver   = Clarabel.Solver(settings)
            Clarabel.setup!(solver, P, q, A, b, cones)
            global sol = Clarabel.solve!(solver)
        """)

        self.solution = jl.sol
        if jl.sol.status == jl.Clarabel.ALMOST_SOLVED:
            warnings.warn(
                "CuClarabel terminated ALMOST_SOLVED; the SDP bound holds only to its "
                "reduced tolerances.",
                stacklevel=2,
            )
        elif jl.sol.status != jl.Clarabel.SOLVED:
            raise RuntimeError(
                f"CuClarabel terminated with {jl.sol.status}; no accepted SDP bound."
            )
        bound = -float(jl.sol.obj_val_dual)
        if not np.isfinite(bound):
            raise RuntimeError("CuClarabel returned a nonfinite SDP dual objective.")
        return bound

__init__

__init__(solver_params=None)

Parameters:

Name	Type	Description	Default
solver_params	dict	Dictionary of solver-specific parameters to configure.	None

Raises:

Type	Description
ValueError	If CUDA dependencies (cupy) are not available.

Source code in fomtools/solvers/cuclarabel.py

def __init__(self, solver_params=None):
    """
    Initialize the GPU-accelerated Clarabel solver.

    Parameters
    ----------
    solver_params : dict, optional
        Dictionary of solver-specific parameters to configure.

    Raises
    ------
    ValueError
        If CUDA dependencies (cupy) are not available.
    """
    _check_cuda_available()
    super(CuClarabelSolver, self).__init__(solver_params=solver_params)

solve

solve(**solve_kwargs)

Solve the SDP on GPU and return the objective value.

Transfers the problem data to GPU memory using CuPy and invokes Julia's CUDA-enabled Clarabel solver via juliacall. Uses cuDSS for the direct linear system solves.

Parameters:

Name	Type	Description	Default
**solve_kwargs		Accepted for API parity, but native Clarabel settings are not yet wired into the Julia GPU path. Passing any option raises ValueError.	{}

Returns:

Type	Description
float	The optimal objective value of the SDP maximization problem.

Notes

Requires CUDA-capable GPU and properly configured Julia environment with CUDA.jl and Clarabel.jl packages.

Source code in fomtools/solvers/cuclarabel.py

def solve(self, **solve_kwargs):
    """
    Solve the SDP on GPU and return the objective value.

    Transfers the problem data to GPU memory using CuPy and invokes
    Julia's CUDA-enabled Clarabel solver via juliacall. Uses cuDSS
    for the direct linear system solves.

    Parameters
    ----------
    **solve_kwargs
        Accepted for API parity, but native Clarabel settings are
        not yet wired into the Julia GPU path. Passing any option
        raises ``ValueError``.

    Returns
    -------
    float
        The optimal objective value of the SDP maximization problem.

    Notes
    -----
    Requires CUDA-capable GPU and properly configured Julia environment
    with CUDA.jl and Clarabel.jl packages.
    """
    merged = {**self.solver_params, **solve_kwargs}
    if merged:
        raise ValueError(
            "cuClarabel does not support custom solver options yet; "
            f"got {merged}. (Construction-time params are rejected up front "
            "by the SDP canonicalizer.)"
        )
    q = self.q
    b_psd = [n * (n + 1) // 2 for n in self.cone_dims["psd"]]
    b = np.concatenate([self.beq, self.bineq, np.zeros(sum(b_psd))])
    A = spa.vstack([self.Aeq, self.Aineq, self.M])

    ndim = q.shape[0]

    import cupy
    from cupyx.scipy.sparse import csr_matrix as cucsr_matrix
    from juliacall import Main as jl

    jl.seval("using CUDA")
    jl.seval("using Clarabel")
    jl.seval("using LinearAlgebra, SparseArrays, CUDA.CUSPARSE")

    # 2. Move Data to GPU
    qgpu = cupy.array(q)
    Agpu = cucsr_matrix(A)
    bgpu = cupy.array(b)

    jl.seval("global q; global A; global b")

    # 5. Transfer Pointers to Julia
    # Note: We initialize P as an empty GPU matrix since it is Zero
    jl.seval(f"""
        P = CuSparseMatrixCSR(sparse(Float64[], Float64[], Float64[], {ndim}, {ndim}))
    """)

    jl.q = jl.Clarabel.cupy_to_cuvector(jl.Float64, int(qgpu.data.ptr), qgpu.size)

    jl.A = jl.Clarabel.cupy_to_cucsrmat(
        jl.Float64,
        int(Agpu.data.data.ptr),
        int(Agpu.indices.data.ptr),
        int(Agpu.indptr.data.ptr),
        *Agpu.shape,
        Agpu.nnz,
    )
    jl.b = jl.Clarabel.cupy_to_cuvector(jl.Float64, int(bgpu.data.ptr), bgpu.size)

    # 6. Setup Cones
    jl.seval("cones = Clarabel.SupportedCone[]")
    jl.push_b(jl.cones, jl.Clarabel.ZeroConeT(self.beq.shape[0]))
    jl.push_b(jl.cones, jl.Clarabel.NonnegativeConeT(self.bineq.shape[0]))
    for n in self.cone_dims["psd"]:
        jl.push_b(jl.cones, jl.Clarabel.PSDTriangleConeT(n))

    # 7. Solve
    # The 'MethodError' will vanish because the Extension now provides the GPU signature.
    # settings = Clarabel.Settings(direct_solve_method = :cudss, verbose = true)
    jl.seval("""
        settings = Clarabel.Settings(direct_solve_method = :cudss, verbose = false)
        solver   = Clarabel.Solver(settings)
        Clarabel.setup!(solver, P, q, A, b, cones)
        global sol = Clarabel.solve!(solver)
    """)

    self.solution = jl.sol
    if jl.sol.status == jl.Clarabel.ALMOST_SOLVED:
        warnings.warn(
            "CuClarabel terminated ALMOST_SOLVED; the SDP bound holds only to its "
            "reduced tolerances.",
            stacklevel=2,
        )
    elif jl.sol.status != jl.Clarabel.SOLVED:
        raise RuntimeError(
            f"CuClarabel terminated with {jl.sol.status}; no accepted SDP bound."
        )
    bound = -float(jl.sol.obj_val_dual)
    if not np.isfinite(bound):
        raise RuntimeError("CuClarabel returned a nonfinite SDP dual objective.")
    return bound