From 4e1798f573ff72dd0f34c9f33b4bc954dbfb945d Mon Sep 17 00:00:00 2001 From: Axel von Kamp Date: Thu, 6 Aug 2026 13:39:02 +0200 Subject: [PATCH] optlang interface for the HiGHS solver supports LP and QP, but QP performance is as good as OSQP --- src/optlang/__init__.py | 8 +- src/optlang/highs_interface.py | 1384 +++++++++++++++++++++ src/optlang/tests/test_highs_interface.py | 536 ++++++++ tox.ini | 2 +- 4 files changed, 1928 insertions(+), 2 deletions(-) create mode 100644 src/optlang/highs_interface.py create mode 100644 src/optlang/tests/test_highs_interface.py diff --git a/src/optlang/__init__.py b/src/optlang/__init__.py index 814ced7..2ff9b19 100644 --- a/src/optlang/__init__.py +++ b/src/optlang/__init__.py @@ -36,6 +36,12 @@ except Exception: log.error('GLPK is available but could not load with error:\n ' + str(traceback.format_exc()).strip().replace('\n','\n ')) +if available_solvers['HIGHS']: + try: + from optlang import highs_interface + except Exception: + log.error('HiGHS is available but could not load with error:\n ' + str(traceback.format_exc()).strip().replace('\n','\n ')) + if available_solvers['CPLEX']: try: from optlang import cplex_interface @@ -69,7 +75,7 @@ # Go through and find the best solver that loaded. Load that one as the default -for engine_str in ['cplex_interface', 'gurobi_interface', 'glpk_interface', +for engine_str in ['cplex_interface', 'gurobi_interface', 'glpk_interface', 'highs_interface', 'hybrid_interface', 'scipy_interface', 'coinor_cbc_interface']: # Must check globals since not all interface variables will be defined if engine_str in globals(): diff --git a/src/optlang/highs_interface.py b/src/optlang/highs_interface.py new file mode 100644 index 0000000..358f8e6 --- /dev/null +++ b/src/optlang/highs_interface.py @@ -0,0 +1,1384 @@ +""" +optlang interface for the HiGHS solver, using the highspy Python bindings. + +Supports LP and QP (Quadratic Programming). +- Continuous variables only (no MILP). +- Linear constraints. +- Linear or Quadratic objectives. + +Install with: + pip install highspy +""" + +import logging +from collections import defaultdict + +import numpy as np + +try: + import highspy +except ImportError: + raise ImportError("The highs_interface requires highspy: pip install highspy") + +from optlang import interface +from optlang import symbolics +from optlang.expression_parsing import parse_optimization_expression +from optlang.exceptions import ContainerAlreadyContains + +log = logging.getLogger(__name__) + + +_HIGHS_STATUS_TO_STATUS = { + "Optimal": interface.OPTIMAL, + "Infeasible": interface.INFEASIBLE, + "Unbounded": interface.UNBOUNDED, + "Unbounded or infeasible": interface.INFEASIBLE_OR_UNBOUNDED, + "Infeasible or unbounded": interface.INFEASIBLE_OR_UNBOUNDED, + "Time limit reached": interface.TIME_LIMIT, + "Iteration limit reached": interface.ITERATION_LIMIT, + "Solve error": interface.UNDEFINED, + "Not Set": interface.UNDEFINED, +} + +_STATUSES_WITH_USABLE_SOLUTION = frozenset([ + interface.OPTIMAL, interface.SUBOPTIMAL, interface.TIME_LIMIT, interface.ITERATION_LIMIT, + interface.INFEASIBLE # needed for cobrapy, see cobra/util/solver.has_primals +]) + + +def _linear_expression_to_dict(expression): + """Turn a linear sympy/symengine expression into ({var_name: coeff}, constant).""" + if getattr(expression, "is_Add", False): + terms = expression.args + else: + terms = (expression,) + + coeffs = defaultdict(float) + constant = 0.0 + + for term in terms: + if term.is_Number: + constant += float(term) + elif term.is_Symbol: + coeffs[term.name] += 1.0 + elif term.is_Mul and len(term.args) == 2: + coeff, sym = term.args + if coeff.is_Number and sym.is_Symbol: + coeffs[sym.name] += float(coeff) + else: + return _linear_expression_to_dict_fallback(expression) + else: + return _linear_expression_to_dict_fallback(expression) + + return dict(coeffs), constant + +def _linear_expression_to_dict_fallback(expression): + offset, linear_coeffs, _ = parse_optimization_expression(None, linear=True, expression=expression) + coeffs = {var.name: float(coeff) for var, coeff in linear_coeffs.items()} + return coeffs, float(offset) + +def _separate_linear_and_quadratic_from_expr(expression): + """ + Separates a raw SymPy expression into linear and quadratic coefficients. + """ + # if isinstance(expression, (int, float)): + # expression = symbolics.sympify(expression) + # else: + expression = expression.expand() + + # 1. Extract Quadratic Terms and get the remaining (linear) expression + quadratic_coeffs, remaining_expr = _get_quadratic_terms_from_expr(expression) + + # 2. Parse the remaining expression using the dedicated linear parser + linear_coeffs, constant = _linear_expression_to_dict(remaining_expr) + + return linear_coeffs, quadratic_coeffs, constant + +def _get_quadratic_terms_from_expr(expression): + """Extract quadratic terms from a raw expression. + + Returns: + quadratic_coeffs (dict): {(var1_name, var2_name): coeff} + remaining_expression (sympy.Expr): The expression with quadratic terms removed. + """ + quadratic_coeffs = defaultdict(float) + + # # Ensure we are working with a SymPy object + # if isinstance(expression, (int, float)): + # expression = symbolics.sympify(expression) + # else: + # expression = expression.expand() + + # We will build the remaining expression by subtracting quadratic terms + # However, subtracting SymPy objects can be tricky with types. + # Instead, we can iterate through the terms and pick out what's NOT quadratic. + + remaining_terms = [] + + def process_term(term): + coeff = 1.0 + factors = [] + + if term.is_Number: + remaining_terms.append(term) + return True + + if term.is_Mul: + args = term.args + else: + args = [term] + + for arg in args: + if arg.is_Number: + coeff *= float(arg) + elif arg.is_Symbol: + factors.append(arg.name) + elif arg.is_Pow: + base, exp = arg.args + if exp.is_Number and float(exp) == 2.0 and base.is_Symbol: + factors.append(base.name) + factors.append(base.name) + else: + # Non-quadratic power, keep it in remaining + remaining_terms.append(term) + return False + else: + # Complex term, keep it + remaining_terms.append(term) + return False + + if len(factors) == 2: + # Quadratic term x*y + v1, v2 = sorted(factors) + quadratic_coeffs[(v1, v2)] += coeff + return True + elif len(factors) == 0: + # This was just a number, handled by is_Number check above + return True + elif len(factors) == 1: + # Linear term + remaining_terms.append(term) + return True + else: + # Degree > 2 + raise ValueError("Degree > 2 detected.") + + if getattr(expression, "is_Add", False): + for term in expression.args: + process_term(term) + else: + process_term(expression) + + # Reconstruct the remaining expression from the kept terms + remaining_expression = symbolics.add(remaining_terms) if remaining_terms else symbolics.sympify(0) + + return dict(quadratic_coeffs), remaining_expression + + +class Variable(interface.Variable): + def __init__(self, name, *args, **kwargs): + if kwargs.get("type", "continuous") != "continuous": + raise ValueError( + "This HiGHS interface only supports " + "continuous variables - LP/QP only, no MIP." + ) + super(Variable, self).__init__(name, **kwargs) + # The column index of this variable in the HiGHS problem it currently + # belongs to (None if it isn't attached to a model). New columns are + # always appended at the end, so Model._add_variables can just read + # off problem.getNumCol() - 1 to set this; on removal, the remaining + # variables' indices are shifted down to match HiGHS's own reindexing + # (see Model._remove_variables). This makes a separate + # name -> index map on the Model unnecessary. + self._solver_index = None + + @interface.Variable.lb.setter + def lb(self, value): + interface.Variable.lb.fset(self, value) + if self.problem is not None: + self.problem._highs_set_col_bounds(self) + + @interface.Variable.ub.setter + def ub(self, value): + interface.Variable.ub.fset(self, value) + if self.problem is not None: + self.problem._highs_set_col_bounds(self) + + @interface.Variable.name.setter + def name(self, value): + # getattr(self, 'problem', None), not self.problem: sympy's + # Symbol.__new__ (via __xnew__) assigns .name during construction, + # before Variable.__init__ has set self.problem at all - a direct + # attribute access would raise AttributeError at that point. + interface.Variable.name.fset(self, value) + if getattr(self, 'problem', None) is not None: + self.problem._highs_set_col_name(self) + + def set_bounds(self, lb, ub): + super(Variable, self).set_bounds(lb, ub) + if self.problem is not None: + self.problem._highs_set_col_bounds(self) + + @property + def primal(self): + if self.problem is None: + return None + return self.problem._variable_primal(self) + + @property + def dual(self): + if self.problem is None: + return None + return self.problem._variable_dual(self) + + +class Constraint(interface.Constraint): + """ + Unlike the earlier version of this interface, this Constraint does not keep + its own copy of the linear coefficients around. Once the constraint is part + of a model, HiGHS's own row data is the single source of truth - writes go + straight through to the solver (mirroring the GLPK interface's design), so + there's no local coefficient cache that can silently drift from it. + + .expression is the one exception: rebuilding it from HiGHS's row data on + every access is expensive, so it's cached in `_expression` and only + rebuilt when `_expression_expired` is set. Anything that changes what a + constraint's row actually contains - a direct coefficient write, or a + referenced variable disappearing out from under it - is responsible for + setting that flag. + + `_initial_coeffs` is only a bootstrap value: it is parsed once at + construction time so `_add_constraints` has something to build the HiGHS + row from, and is never consulted again once that row exists. + + `_constant` remains local state (HiGHS rows don't carry a constant offset + themselves, so it has to be folded into the row bounds), but a bare + constant can never reference a variable, so it can never go stale the way + a coefficient cache can. + """ + + _INDICATOR_CONSTRAINT_SUPPORT = False + + def __init__(self, expression, sloppy=False, *args, **kwargs): + if isinstance(expression, (int, float)): + sloppy = True + expression = symbolics.Real(expression) + + self._expression_expired = False + super(Constraint, self).__init__(expression, sloppy=sloppy, *args, **kwargs) + + if not sloppy and not self.is_Linear: + raise ValueError( + "The HiGHS interface only supports linear constraints. " + "%s is not linear." % self + ) + + self._initial_coeffs, self._constant = _linear_expression_to_dict(self.expression) + # The row index of this constraint in the HiGHS problem it currently + # belongs to (None if it isn't attached to a model, or not yet added + # to the model's HiGHS instance). See Variable._solver_index for why + # this replaces a name -> index map kept on the Model. + self._solver_index = None + + @interface.Constraint.lb.setter + def lb(self, value): + interface.Constraint.lb.fset(self, value) + if self.problem is not None: + self.problem._highs_set_row_bounds(self) + + @interface.Constraint.ub.setter + def ub(self, value): + interface.Constraint.ub.fset(self, value) + if self.problem is not None: + self.problem._highs_set_row_bounds(self) + + @interface.Constraint.name.setter + def name(self, value): + interface.Constraint.name.fset(self, value) + if getattr(self, 'problem', None) is not None: + self.problem._highs_set_row_name(self) + + @property + def primal(self): + if self.problem is None: + return None + return self.problem._constraint_primal(self) + + @property + def dual(self): + if self.problem is None: + return None + return self.problem._constraint_dual(self) + + def _get_expression(self): + """Reconstruct the expression by reading HiGHS's live row data if necessary.""" + if not self._expression_expired: + return self._expression + elif self.problem is not None: # constraint may be pending removal + self.problem.update() + + if self.problem is None: + raise ValueError("Constraint expression cannot be constructed.") + + row_index = self._solver_index + if row_index is None: + # Constructed, but not yet added to the model's HiGHS instance - + # report the bootstrap coefficients that _add_constraints will use. + terms = [] + for var_name, coeff in self._initial_coeffs.items(): + var = self.problem.variables.get(var_name) + if var is not None: + terms.append(coeff * var) + else: + _, idx, val = self.problem.problem.getRowEntries(row_index) + terms = [] + for i, v in zip(idx, val): + var_name = self.problem.problem.variableName(int(i)) + var = self.problem.variables.get(var_name) + if var is not None: + terms.append(float(v) * var) + + if self._constant != 0: + # symbolics.add() uses SymPy's low-level Add._from_args(), which + # (unlike Add(*args)) does not auto-sympify its arguments - a bare + # Python float here breaks it (AttributeError: 'float' object has + # no attribute 'is_commutative'), so sympify explicitly. + terms.append(symbolics.Real(self._constant)) + + self._expression = symbolics.add(terms) if terms else symbolics.sympify(0) + self._expression_expired = False + return self._expression + + def set_linear_coefficients(self, coefficients, sloppy=False): + """ + Sets the linear coefficients of the constraint directly, bypassing SymPy + parsing. Writes go straight through to HiGHS's row - there's no local + cache to keep in sync, so this can never drift from the solver's state. + + Args: + coefficients (dict): A dictionary of {var: coeff}. + sloppy (bool): If True, skip re-resolving each variable through + the model (self.problem.variables.get(var.name)) and pass + `coefficients` straight through to HiGHS. Only safe when the + variables in `coefficients` are already known to belong to + this constraint's model (e.g. they came from model.variables + in the first place) - saves a dict lookup per variable, which + matters when this is called in a hot loop. + """ + if self.problem is None: + raise Exception("Can't change coefficients if constraint is not associated with a model.") + + if self._solver_index is None: + # Not yet added to the model - update the bootstrap coefficients + # that will be used to build the row. + for var, coeff in coefficients.items(): + self._initial_coeffs[var.name] = coeff + return + + if sloppy: + self.problem._highs_set_coefficients(self, coefficients) + self._expression_expired = True + return + + model_coeffs = {} + for var, coeff in coefficients.items(): + var_obj = self.problem.variables.get(var.name) + if var_obj is not None: + model_coeffs[var_obj] = coeff + + self.problem._highs_set_coefficients(self, model_coeffs) + self._expression_expired = True + + def get_linear_coefficients(self, variables): + """ + Get coefficients of linear terms in constraint. + + Like everything else on this class, this reads straight from HiGHS's + row data when the constraint is attached to a model (via the same + helper _get_expression uses), rather than from any local cache. + + Parameters + ---------- + variables : iterable + An iterable of Variable objects + + Returns + ------- + Coefficients : dict + {var1: coefficient, var2: coefficient ...}. Variables with no + entry in the constraint get a coefficient of 0. + """ + if self.problem is None: + raise Exception( + "Can't get coefficients from solver if constraint is not associated with a model.") + coeffs, _, _ = self.problem._constraint_to_coeffs(self) + return {variable: coeffs.get(variable.name, 0) for variable in variables} + + +class Objective(interface.Objective): + """ + The linear part of the objective follows the same ground-truth design as + Constraint: once attached to a model, HiGHS's own cost vector (readable via + Highs.getObjective()) is authoritative, and coefficient writes go straight + through to Highs.changeColCost. There's no local linear-coefficient cache + left to go stale. + + The quadratic part (_quadratic_coeffs) is still cached locally. highspy's + high-level wrapper (highs.py) has no Hessian read-back equivalent to + getObjective()/getRowEntries() - passHessian is write-only from what's + exposed there. Until that's confirmed one way or another against the + compiled _core module, the Hessian remains the one piece still built from + a local cache rather than queried live, and _remove_variables still has to + prune it explicitly when a variable disappears. + """ + + def __init__(self, expression, sloppy=False, *args, **kwargs): + if isinstance(expression, (int, float)): + expression = symbolics.Real(expression) + sloppy = True + + kwargs['sloppy'] = sloppy + super(Objective, self).__init__(expression, *args, **kwargs) + + if not sloppy: + if not (self.is_Linear or self.is_Quadratic): + raise ValueError("The HiGHS interface only supports linear or quadratic objectives.") + + # Bootstrap values only: _initial_linear_coeffs seeds HiGHS's cost + # vector the first time this Objective is attached to a model (see + # Model._update_linear_objective) and is never consulted again after + # that. _quadratic_coeffs remains the ongoing source of truth for the + # Hessian (see class docstring). + self._initial_linear_coeffs, self._quadratic_coeffs, self._constant = \ + _separate_linear_and_quadratic_from_expr(self.expression) + + def _get_expression(self): + """Reconstruct the expression, reading the linear part live from HiGHS.""" + if self.problem is None: + return self._expression + + obj_expr, _ = self.problem.problem.getObjective() + linear_coeffs = {} + for i, v in zip(obj_expr.idxs, obj_expr.vals): + var_name = self.problem.problem.variableName(int(i)) + linear_coeffs[var_name] = float(v) + + terms = [] + + for (v1_name, v2_name), coeff in self._quadratic_coeffs.items(): + v1 = self.problem.variables.get(v1_name) + v2 = self.problem.variables.get(v2_name) + if v1 is not None and v2 is not None: + terms.append(coeff * v1 * v2) + + for var_name, coeff in linear_coeffs.items(): + var = self.problem.variables.get(var_name) + if var is not None: + terms.append(coeff * var) + + if self._constant != 0: + # See Constraint._get_expression - symbolics.add() needs an + # already-sympified argument, not a bare Python float. + terms.append(symbolics.sympify(self._constant)) + + return symbolics.add(terms) if terms else symbolics.sympify(0) + + @interface.Objective.expression.getter + def expression(self): + return self._get_expression() + + @interface.Objective.direction.setter + def direction(self, value): + interface.Objective.direction.fset(self, value) + if self.problem is not None: + self.problem._highs_set_objective_sense(value) + + @property + def value(self): + if self.problem is None: + return None + return self.problem._get_objective_value() + + def set_linear_coefficients(self, coefficients): + self.set_coefficients(linear=coefficients) + + def set_coefficients(self, linear=None, quadratic=None, constant=None): + """ + Sets the objective coefficients directly, bypassing SymPy parsing. + + Args: + linear (dict, optional): Dictionary of {Variable: coeff} for linear terms. + quadratic (dict, optional): Dictionary of {(Variable, Variable): coeff} + for quadratic terms. + constant (float, optional): Constant term of the objective. + """ + if self.problem is None: + raise Exception("Can't change coefficients if objective is not associated with a model.") + + # Linear coefficients go straight through to HiGHS's cost vector - + # only the named variables are touched, so this is naturally a merge, + # not a replace, with no local bookkeeping required. (This matters in + # practice: callers like cobrapy build up an objective by calling this + # incrementally, one or two variables at a time, across many calls, + # and expect earlier calls' coefficients to survive.) + if linear is not None: + self.problem._highs_set_objective_coefficients(linear) + + # Quadratic coefficients are different: _quadratic_coeffs is still a + # local cache standing in for the whole Hessian (see class + # docstring), and passing quadratic={} is the only way to express + # "there are now fewer quadratic terms than before" - including zero + # of them, e.g. when moving from a QP back to a plain LP. So unlike + # linear, an explicitly-given quadratic dict *replaces* the cache + # outright rather than merging into it. + if quadratic is not None: + self._quadratic_coeffs = { + (v1.name, v2.name): coeff for (v1, v2), coeff in quadratic.items() + } + self.problem._update_quadratic_objective() + + if constant is not None: + self._constant = float(constant) + self.problem.problem.changeObjectiveOffset(float(constant)) + + def get_linear_coefficients(self, variables): + """ + Get coefficients of linear terms in objective. + + Reads the linear part straight from HiGHS's live cost vector when + attached to a model, the same way _get_expression does. The + quadratic part is intentionally not reflected here since this only + reports linear coefficients (see set_linear_coefficients/ + set_coefficients for the quadratic equivalent). + + Parameters + ---------- + variables : iterable + An iterable of Variable objects + + Returns + ------- + Coefficients : dict + {var1: coefficient, var2: coefficient ...}. Variables with no + linear entry in the objective get a coefficient of 0. + """ + if self.problem is None: + raise Exception( + "Can't get coefficients from solver if objective is not associated with a model.") + obj_expr, _ = self.problem.problem.getObjective() + linear_coeffs = { + self.problem.problem.variableName(int(i)): float(v) + for i, v in zip(obj_expr.idxs, obj_expr.vals) + } + return {variable: linear_coeffs.get(variable.name, 0) for variable in variables} + + +class Tolerances(object): + def __init__(self, configuration): + self._configuration = configuration + self._feasibility = 1e-7 + self._optimality = 1e-7 + # HiGHS always solves QPs with its interior-point method, which has its + # own convergence tolerance distinct from the simplex-oriented + # primal/dual feasibility tolerances above. Exposed separately since + # loosening it is often the fix for QPs that fail to converge for + # purely numerical/scaling reasons. + self._ipm_optimality = 1e-8 + + @property + def feasibility(self): + return self._feasibility + + @feasibility.setter + def feasibility(self, value): + self._feasibility = value + self._configuration._apply_tolerances() + + @property + def optimality(self): + return self._optimality + + @optimality.setter + def optimality(self, value): + self._optimality = value + self._configuration._apply_tolerances() + + @property + def ipm_optimality(self): + return self._ipm_optimality + + @ipm_optimality.setter + def ipm_optimality(self, value): + self._ipm_optimality = value + self._configuration._apply_tolerances() + + def to_dict(self): + return { + "feasibility": self.feasibility, + "optimality": self.optimality, + "ipm_optimality": self.ipm_optimality, + } + + +class Configuration(interface.MathematicalProgrammingConfiguration): + def __init__(self, verbosity=0, timeout=None, presolve=True, + problem=None, *args, **kwargs): + self.problem = problem + self._verbosity = verbosity + self._timeout = timeout + self._presolve = presolve + self._tolerances = Tolerances(self) + super(Configuration, self).__init__(*args, **kwargs) + self.problem = problem + if self.problem is not None: + self._apply_all() + + def _apply_all(self): + h = self.problem.problem + if h is None: + return + h.setOptionValue("output_flag", self._verbosity != 0) + if self._timeout is not None: + h.setOptionValue("time_limit", float(self._timeout)) + h.setOptionValue("presolve", "on" if self._presolve else "off") + self._apply_tolerances() + + def _apply_tolerances(self): + if self.problem is not None and self.problem.problem is not None: + h = self.problem.problem + h.setOptionValue("primal_feasibility_tolerance", self._tolerances.feasibility) + h.setOptionValue("dual_feasibility_tolerance", self._tolerances.feasibility) + h.setOptionValue("ipm_optimality_tolerance", self._tolerances.ipm_optimality) + + @property + def tolerances(self): + return self._tolerances + + @tolerances.setter + def tolerances(self, value): + if isinstance(value, Tolerances): + self._tolerances = value + + @property + def verbosity(self): + return self._verbosity + + @verbosity.setter + def verbosity(self, value): + self._verbosity = value + if self.problem is not None and self.problem.problem is not None: + self.problem.problem.setOptionValue("output_flag", value != 0) + + @property + def timeout(self): + return self._timeout + + @timeout.setter + def timeout(self, value): + self._timeout = value + if self.problem is not None and self.problem.problem is not None: + if value is not None: + self.problem.problem.setOptionValue("time_limit", float(value)) + + @property + def presolve(self): + return self._presolve + + @presolve.setter + def presolve(self, value): + self._presolve = value + if self.problem is not None and self.problem.problem is not None: + self.problem.problem.setOptionValue("presolve", "on" if value else "off") + + +class _HighsBackedContainer(object): + """ + Shared base for the model.variables / model.constraints containers. Their + ordering and indexing are read directly from the live HiGHS instance - + the same "ask the solver" principle used throughout this interface for + Constraint/Objective coefficients - instead of a private _object_list / + _indices that could drift out of sync with the solver's own column/row + order (e.g. after HiGHS shifts indices down on _remove_variables's + deleteVars call). + + A name -> object dict is still kept, since HiGHS itself has no notion of + the optlang Variable/Constraint wrapper objects, only of names and + column/row numbers - "how many are there" and "what order are they in" + come from HiGHS, "which Python object goes with this name" comes from + this dict. + + This deliberately does not subclass optlang.container.Container. That + class carries its own private _object_list/_indices bookkeeping (exactly + what this class exists to avoid), plus API surface this interface never + touches: __setitem__, fromkeys, and the py2-era iterkeys/itervalues/ + iteritems/has_key aliases. Implementing just the subset actually used by + optlang.interface.Model and this file - append/extend, __contains__/ + __iter__/__getitem__/__delitem__/__len__, get/keys/values/items, and + update_key (needed by Variable.name's and Constraint.name's setters when + an attached object is renamed) - keeps this class's behavior easy to + verify against what's really exercised. + + Subclasses only need to implement _num() and _name_at(idx). + """ + + def __init__(self, model, iterable=()): + self._model = model + self._dict = {} + for item in iterable: + self.append(item) + + def _num(self): + raise NotImplementedError("Subclasses must implement _num().") + + def _name_at(self, idx): + raise NotImplementedError("Subclasses must implement _name_at().") + + def _live_names(self): + return [self._name_at(i) for i in range(self._num())] + + def __len__(self): + return self._num() + + def __contains__(self, item): + if isinstance(item, str): + return item in self._dict + name = getattr(item, "name", None) + return name is not None and name in self._dict and item is self._dict[name] + + def __iter__(self): + for name in self._live_names(): + obj = self._dict.get(name) + if obj is not None: + yield obj + + def __getitem__(self, item): + if isinstance(item, slice): + names = self._live_names()[item] + return [self._dict[name] for name in names if name in self._dict] + if isinstance(item, int): + names = self._live_names() + return self._dict[names[item]] + return self._dict[item] # by name + + def __delitem__(self, key): + name = key if isinstance(key, str) else self[key].name + del self._dict[name] + + def keys(self): + return list(self._live_names()) + + def values(self): + return list(self) + + def items(self): + return [(name, self._dict[name]) for name in self._live_names() if name in self._dict] + + def get(self, key, default=None): + try: + return self[key] + except (KeyError, IndexError): + return default + + def append(self, value): + if not hasattr(value, "name"): + raise AttributeError( + 'Object %s does not have a "name" attribute and cannot be stored.' % value) + name = value.name + if name in self._dict: + raise ContainerAlreadyContains( + "Container '%s' already contains an object with name '%s'." % (self, value.name)) + self._dict[name] = value + + def extend(self, values): + for value in values: + self.append(value) + + def update_key(self, key): + """Re-files an already-stored item under its current .name. + + Called by Variable.name's and Constraint.name's setters (in + optlang.interface) right after an attached object's name changes, so + this dict's key stays in sync with the object it points at. + """ + item = self._dict.get(key) + if item is None: + return + name = item.name + if key != name: + self._dict[name] = item + del self._dict[key] + + def __getstate__(self): + return list(self) + + def __setstate__(self, obj_list): + # Unlike optlang.container.Container, this can't just replay + # self.__init__(obj_list): a live HiGHS-backed container needs its + # owning Model, which isn't part of the pickled state. + # Reconstructing the Model rebuilds this instead. + raise NotImplementedError( + "%s is backed by a live HiGHS instance and cannot be unpickled on " + "its own - restore the owning Model instead." % type(self).__name__ + ) + + def __getattr__(self, name): + try: + return self[name] + except KeyError: + raise AttributeError("'%s' object has no attribute %s" % (self, name)) + + def __dir__(self): + attributes = list(self.__class__.__dict__.keys()) + attributes.extend(self._dict.keys()) + return attributes + + +class HighsVariablesContainer(_HighsBackedContainer): + """model.variables, backed by HiGHS's live column order/names.""" + + def _num(self): + return self._model.problem.getNumCol() + + def _name_at(self, idx): + return self._model.problem.variableName(idx) + + +class HighsConstraintsContainer(_HighsBackedContainer): + """model.constraints, backed by HiGHS's live row order/names.""" + + def _num(self): + return self._model.problem.getNumRow() + + def _name_at(self, idx): + _, name = self._model.problem.getRowName(idx) + return name + + +class Model(interface.Model): + def __init__(self, problem=None, *args, **kwargs): + self.problem = None + self._has_solution = False + self._solution_col_value = None + self._solution_col_dual = None + self._solution_row_value = None + self._solution_row_dual = None + self._objective_value = None + self._last_objective_was_quadratic = False + super(Model, self).__init__(problem=problem, *args, **kwargs) + + def _initialize_problem(self): + self.problem: highspy.Highs = highspy.Highs() + self.problem.setOptionValue("output_flag", False) + # Replace the base class's plain Containers (created before this + # point, so still empty here) with ones backed by this Highs + # instance's live column/row order. + self._variables = HighsVariablesContainer(self) + self._constraints = HighsConstraintsContainer(self) + + def _initialize_model_from_problem(self, problem: highspy.Highs): + """ + Copies an already-populated highspy.Highs instance, reconstructing + optlang Variable/Constraint/Objective wrappers for its current + columns, rows and cost vector. + """ + if not isinstance(problem, highspy.Highs): + raise TypeError( + "Provided problem must be a highspy.Highs instance, not %s." % type(problem)) + + # self.problem = problem + # self._variables = HighsVariablesContainer(self) + # self._constraints = HighsConstraintsContainer(self) + + # lp = problem.getLp() + + self._initialize_problem() # create a new highspy.Highs problem + self.problem.passModel(problem.getModel()) # copy the source problem + lp = self.problem.getLp() + problem = self.problem + + inf = highspy.kHighsInf + + # Variables: one wrapper per existing column, bounds and index taken + # straight from the live problem. + for i in range(problem.getNumCol()): + name = problem.variableName(i) + lb = float(lp.col_lower_[i]) + ub = float(lp.col_upper_[i]) + variable = Variable( + name, + lb=None if lb <= -inf else lb, + ub=None if ub >= inf else ub, + ) + # Raw assignment, not variable.problem = self: mirrors + # _add_variables/the rest of this file, and (unlike Constraint's + # or Objective's `problem`) Variable.problem is a plain + # attribute with no side effects, so either would work - this + # just keeps the pattern identical across all three. + variable._solver_index = i + variable.problem = self + self._variables.append(variable) + self._variables_to_constraints_mapping[name] = set() + + # Constraints: one wrapper per existing row. Built the same way + # Constraint.__init__ always is elsewhere in this file - constructed + # unattached, with _solver_index and _problem only patched in + # afterwards. Constraint._get_expression touches self._solver_index, + # so constructing it directly with problem=self here (before that + # attribute exists) would blow up; this order avoids that. + for i in range(problem.getNumRow()): + _, name = problem.getRowName(i) + _, row_lb, row_ub, _ = problem.getRow(i) + _, idx, val = problem.getRowEntries(i) + # terms = [float(v) * self._variables[problem.variableName(int(j))] for j, v in zip(idx, val)] + # expression = symbolics.add(terms) if terms else symbolics.sympify(0) + # constraint = Constraint( + constraint = Constraint( + 0, # do not set up explicit expression here, can be constructed on demand + lb=None if row_lb <= -inf else float(row_lb), + ub=None if row_ub >= inf else float(row_ub), + name=name, + sloppy=True, + ) + constraint._expression_expired = True + constraint._solver_index = i + constraint._problem = self + self._constraints.append(constraint) + for j in idx: + var_name = problem.variableName(int(j)) + self._variables_to_constraints_mapping[var_name].add(name) + + # Objective: linear part only, see docstring. Same construction + # order rationale as Constraint above - Objective._get_expression + # reads self._quadratic_coeffs, which does not exist until after + # __init__ returns, so this too is built unattached first. + obj_expr, sense = problem.getObjective() + linear_coeffs = {} + terms = [] + for i, v in zip(obj_expr.idxs, obj_expr.vals): + var_name = problem.variableName(int(i)) + linear_coeffs[var_name] = float(v) + terms.append(float(v) * self._variables[var_name]) + constant = float(obj_expr.constant or 0.0) + if constant: + terms.append(symbolics.Real(constant)) + expression = symbolics.add(terms) if terms else symbolics.sympify(0) + direction = "max" if sense == highspy.ObjSense.kMaximize else "min" + + objective = Objective(expression, direction=direction, sloppy=True) + objective._initial_linear_coeffs = linear_coeffs + objective._quadratic_coeffs = {} + objective._constant = constant + objective._problem = self + self._objective = objective + + def _build_hessian_data(self, quadratic_coeffs, num_col): + """ + Helper to build Hessian matrix data (start, index, value) for HiGHS. + Returns (q_start, q_index, q_value) or (None, None, None) if empty. + """ + if not quadratic_coeffs: + return None, None, None + + q_start = np.zeros(num_col + 1, dtype=np.int32) + q_index_list = [] + q_value_list = [] + q_col_entries = [[] for _ in range(num_col)] + + for (v1_name, v2_name), coeff in quadratic_coeffs.items(): + v1 = self._variables.get(v1_name) + v2 = self._variables.get(v2_name) + if v1 is None or v2 is None: + continue + i = v1._solver_index + j = v2._solver_index + + # Lower Triangular (i >= j) + if i < j: + i, j = j, i + + val = coeff + if i == j: + val = 2.0 * coeff + + q_col_entries[j].append((i, val)) + + # Sort and build CSC arrays + pos = 0 + for j in range(num_col): + q_start[j] = pos + q_col_entries[j].sort(key=lambda x: x[0]) + for (row_idx, value) in q_col_entries[j]: + q_index_list.append(row_idx) + q_value_list.append(value) + pos += 1 + q_start[num_col] = pos + + return q_start, np.array(q_index_list, dtype=np.int32), np.array(q_value_list, dtype=np.double) + + def _update_linear_objective(self): + """ + Fully (re)seeds the linear part (costs, offset, sense) of the objective + in HiGHS. This only runs when the whole Objective object is (re)assigned + via the `model.objective = ...` setter, so wiping all costs first and + seeding from the new objective's bootstrap coefficients is intentional - + it's a full replace, not an incremental update. Incremental updates go + through Objective.set_coefficients / _highs_set_objective_coefficients + instead, which only ever touch the specific columns being changed. + """ + if self.problem is None or self.objective is None: + return + + n = self.problem.getNumCol() + if n > 0: + self.problem.changeColsCost(n, np.arange(n, dtype=np.int32), np.zeros(n, dtype=np.double)) + + # Re-derive from the objective's own (unattached-side) expression + # rather than trusting _initial_linear_coeffs/_constant verbatim: + # those are snapshotted once at Objective.__init__ time, but + # +=/-=/*= (OptimizationExpression's __iadd__/__isub__/__imul__) + # mutate _expression directly and have no way to know they should + # keep that snapshot in sync. Re-parsing here, at the one place the + # snapshot is actually consumed, keeps a mutated-before-attachment + # Objective (e.g. `obj = Objective(x); obj += 2 * y; model.objective + # = obj`) working correctly without having to intercept every + # mutating operator individually. + linear_coeffs, _, obj_constant = _separate_linear_and_quadratic_from_expr(self.objective._expression) + self.objective._initial_linear_coeffs = linear_coeffs + self.objective._constant = obj_constant + + for var_name, coeff in linear_coeffs.items(): + var = self._variables.get(var_name) + if var is not None: + self.problem.changeColCost(var._solver_index, float(coeff)) + + self.problem.changeObjectiveOffset(float(obj_constant)) + self._highs_set_objective_sense(self.objective.direction) + + def _update_quadratic_objective(self): + """Updates the quadratic part (Hessian) of the objective in HiGHS.""" + if self.problem is None or self.objective is None: + return + + # Same rationale as the _highs_set_* helpers below: quadratic terms + # may reference variables added but not yet flushed. + self.update() + + quadratic_coeffs = self.objective._quadratic_coeffs + num_col = self.problem.getNumCol() + + q_start, q_index, q_value = self._build_hessian_data(quadratic_coeffs, num_col) + + if q_start is not None: + hessian = highspy.HighsHessian() + hessian.dim_ = num_col + hessian.format_ = highspy.HessianFormat.kTriangular + hessian.start_ = q_start + hessian.index_ = q_index + hessian.value_ = q_value + self.problem.passHessian(hessian) + else: + # Clear the Hessian if no quadratic terms + empty_hessian = highspy.HighsHessian() + empty_hessian.dim_ = num_col + empty_hessian.format_ = highspy.HessianFormat.kTriangular + empty_hessian.start_ = np.zeros(num_col + 1, dtype=np.int32) + empty_hessian.index_ = np.array([], dtype=np.int32) + empty_hessian.value_ = np.array([], dtype=np.double) + self.problem.passHessian(empty_hessian) + + self._last_objective_was_quadratic = len(quadratic_coeffs) > 0 + + def _sync_objective_to_solver(self): + """Sync the current objective's internal state to the HiGHS solver.""" + self._update_linear_objective() + self._update_quadratic_objective() + + @interface.Model.objective.setter + def objective(self, value): + interface.Model.objective.fset(self, value) + if self.problem is not None: + self._sync_objective_to_solver() + + def _add_variables(self, variables): + super(Model, self)._add_variables(variables) + inf = highspy.kHighsInf + for variable in variables: + lb = -inf if variable.lb is None else variable.lb + ub = inf if variable.ub is None else variable.ub + # name=... is required here: every read path (Constraint/Objective + # .expression, _constraint_to_coeffs) maps a + # HiGHS column index back to an optlang variable name via + # variableName(idx), which raises if the column was never named. + self.problem.addVariable(lb, ub, name=variable.name) + # New columns are always appended at the end, so the newly added + # variable's index is simply the last column. + variable._solver_index = self.problem.getNumCol() - 1 + + def _add_constraints(self, constraints, sloppy=False): + super(Model, self)._add_constraints(constraints, sloppy=sloppy) + inf = highspy.kHighsInf + for constraint in constraints: + coeffs, lb, ub = self._constraint_to_coeffs(constraint) + row_lb = -inf if lb is None else lb + row_ub = inf if ub is None else ub + indices = np.array([self._variables[name]._solver_index for name in coeffs.keys()], dtype=np.int32) + values = np.array(list(coeffs.values()), dtype=np.double) + self.problem.addRow(row_lb, row_ub, len(indices), indices, values) + # New rows are always appended at the end, so the newly added + # constraint's index is simply the last row. + row_index = self.problem.getNumRow() - 1 + constraint._solver_index = row_index + # addRow (unlike addVariable) has no name kwarg, so name the row + # in a separate call. Needed by variableName's row counterpart + # (getRowName) used in _constraint_to_coeffs and + # HighsConstraintsContainer. + self.problem.passRowName(row_index, constraint.name) + + def _remove_variables(self, variables): + for variable in variables: + try: + self._variables[variable.name] + except KeyError: + raise LookupError("Variable %s not in solver" % variable.name) + + removed_names = set(variable.name for variable in variables) + removed_indices = sorted(variable._solver_index for variable in variables) + old_num_col = self.problem.getNumCol() + + # Physically delete the columns from HiGHS, not just from our own + # bookkeeping. This is essential under the ground-truth design: + # Constraint/Objective now read their linear coefficients straight + # from HiGHS's live row/cost data rather than from a local cache, so + # if the column itself is left in place, every row that referenced it + # keeps reporting it - even though our _variables container has + # already forgotten about it. Deleting the column here lets HiGHS + # itself cascade the removal through every row and the objective, + # exactly like GLPK's glp_del_cols. + self.problem.deleteVars(len(removed_indices), removed_indices) + + # HiGHS shifts every remaining column's index down to close the gaps + # left by the deleted ones - shift the remaining variables' own + # _solver_index attributes to match. But if every removed column was + # already at the tail (the common case: variables added and then + # removed again within the same `with model:` block, e.g. slack + # variables in core.make_scenario_feasible), there are no gaps to + # close below them and nothing to shift - skip the loop entirely. + # (Sorted, unique indices whose smallest value equals + # old_num_col - len(removed_indices) can only be the top + # len(removed_indices) indices; there's no other way to fit that + # many distinct values into that range.) + if removed_indices[0] != old_num_col - len(removed_indices): + for variable in self._variables: + if variable.name in removed_names: + continue + shift = sum(1 for r in removed_indices if r < variable._solver_index) + variable._solver_index -= shift + + for variable in variables: + # Constraint now caches its expression (_expression_expired) as + # of the sloppy-removal changes below - a constraint referencing + # this variable has a stale cache the moment the variable's + # column disappears from its row, so mark it expired before + # forgetting the mapping that's the only way to find it. (The + # comment this replaced predated that cache and was no longer + # accurate: it used to be true that there was nothing to + # invalidate here, back when Constraint.expression always read + # HiGHS live with no local cache at all.) + for constraint_name in self._variables_to_constraints_mapping.get(variable.name, ()): + constraint = self._constraints.get(constraint_name) + if constraint is not None: + constraint._expression_expired = True + self._variables_to_constraints_mapping.pop(variable.name, None) + variable.problem = None + del self._variables[variable.name] + + # The objective's quadratic part is still cached locally (see + # Objective's docstring - no live Hessian read-back is available yet), + # so it's the one remaining piece that needs explicit pruning: a + # lingering quadratic term referencing a removed variable would break + # _build_hessian_data the next time the objective is rebuilt. + if self.objective is not None and self.objective._quadratic_coeffs: + self.objective._quadratic_coeffs = { + k: v for k, v in self.objective._quadratic_coeffs.items() + if k[0] not in removed_names and k[1] not in removed_names + } + + def _remove_constraints(self, constraints): + for constraint in constraints: + try: + self._constraints[constraint.name] + except KeyError: + raise LookupError("Constraint %s not in solver" % constraint) + + removed_names = set(constraint.name for constraint in constraints) + removed_indices = sorted(constraint._solver_index for constraint in constraints) + old_num_row = self.problem.getNumRow() + + # Mirrors _remove_variables: physically delete the rows from HiGHS + # and let it cascade the reindexing, rather than rebuilding the + # whole LP from scratch. Removing rows doesn't touch any column, so + # unlike variable removal there's no Hessian/objective cache to + # prune here - constraints don't appear in any other object's cache. + self.problem.deleteRows(len(removed_indices), removed_indices) + + # Same end-of-range shortcut as _remove_variables: if every removed + # row was already at the tail, nothing below them moves. + if removed_indices[0] != old_num_row - len(removed_indices): + for constraint in self._constraints: + if constraint.name in removed_names: + continue + shift = sum(1 for r in removed_indices if r < constraint._solver_index) + constraint._solver_index -= shift + + for constraint in constraints: + constraint.problem = None + del self._constraints[constraint.name] + + def _constraint_bounds(self, constraint: Constraint): + lb, ub = constraint.lb, constraint.ub + if lb is not None: + lb = lb - constraint._constant + if ub is not None: + ub = ub - constraint._constant + return lb, ub + + def _constraint_to_coeffs(self, constraint: Constraint): + lb, ub = self._constraint_bounds(constraint) + + row_index = constraint._solver_index + if row_index is not None: + # The row already exists in HiGHS, which is the source of truth + # for its own coefficients - read them back directly rather than + # trusting any local copy (there isn't one anymore anyway). + _, idx, val = self.problem.getRowEntries(row_index) + coeffs = { + self.problem.variableName(int(i)): float(v) + for i, v in zip(idx, val) + } + else: + # Constraint not yet added to this model - use its bootstrap + # coefficients (this is the path _add_constraints takes). + coeffs = dict(constraint._initial_coeffs) + + return coeffs, lb, ub + + def _highs_set_col_bounds(self, variable): + # Flush any pending add()s first: a variable can have its lb/ub + # setter (or set_bounds) called right after being added to a model, + # before update() has run to assign it a _solver_index. Without this, + # changeColBounds would be called with index None. update() is a + # cheap no-op when there is nothing pending, so this costs nothing + # in the common case where the variable was already committed. + self.update() + inf = highspy.kHighsInf + lb = -inf if variable.lb is None else variable.lb + ub = inf if variable.ub is None else variable.ub + self.problem.changeColBounds(variable._solver_index, lb, ub) + + def _highs_set_row_bounds(self, constraint): + # See _highs_set_col_bounds -- same rationale, for constraint rows. + self.update() + inf = highspy.kHighsInf + lb, ub = self._constraint_bounds(constraint) + row_lb = -inf if lb is None else lb + row_ub = inf if ub is None else ub + self.problem.changeRowBounds(constraint._solver_index, row_lb, row_ub) + + def _highs_set_col_name(self, variable): + # Same flush rationale as _highs_set_col_bounds. + self.update() + self.problem.passColName(variable._solver_index, variable.name) + + def _highs_set_row_name(self, constraint): + # Same flush rationale as _highs_set_col_bounds. + self.update() + self.problem.passRowName(constraint._solver_index, constraint.name) + + def _highs_set_coefficients(self, constraint, coefficients): + # Same rationale as _highs_set_col_bounds: coefficients may reference + # a variable that was just add()ed but not yet update()d, so its + # _solver_index would still be None. Constraint.set_linear_coefficients + # already resolves variables through the auto-updating `variables` + # property in its default (non-sloppy) path, making this redundant + # there -- but the sloppy=True path bypasses that resolution + # entirely and calls straight through to here, so this is the one + # place that also has to hold for that path to be safe by default. + self.update() + i = constraint._solver_index + for var, coeff in coefficients.items(): + j = var._solver_index + self.problem.changeCoeff(i, j, float(coeff)) + + def _highs_set_objective_coefficients(self, coefficients): + # Unlike Constraint.set_linear_coefficients' default path, Objective's + # linear coefficient path does not resolve variables through the + # auto-updating `variables` property before reaching here -- so this + # is the one place that needs to force freshness for it. + self.update() + for var, coeff in coefficients.items(): + self.problem.changeColCost(var._solver_index, float(coeff)) + + def _highs_set_objective_sense(self, direction): + self.problem.changeObjectiveSense( + highspy.ObjSense.kMaximize if direction == "max" else highspy.ObjSense.kMinimize + ) + + def _optimize(self): + h = self.problem + h.run() + + status_str = h.modelStatusToString(h.getModelStatus()) + status = _HIGHS_STATUS_TO_STATUS.get(status_str, interface.UNDEFINED) + + self._has_solution = status in _STATUSES_WITH_USABLE_SOLUTION + self._solution_col_value = None + self._solution_col_dual = None + self._solution_row_value = None + self._solution_row_dual = None + self._objective_value = None + + return status + + def _ensure_solution_arrays(self): + if self._solution_col_value is not None: + return + solution = self.problem.getSolution() + self._solution_col_value = np.asarray(list(solution.col_value)) + self._solution_col_dual = np.asarray(list(solution.col_dual)) + self._solution_row_value = np.asarray(list(solution.row_value)) + self._solution_row_dual = np.asarray(list(solution.row_dual)) + + def _variable_primal(self, variable): + if not self._has_solution: + return None + self._ensure_solution_arrays() + return float(self._solution_col_value[variable._solver_index]) + + def _variable_dual(self, variable): + if not self._has_solution: + return None + self._ensure_solution_arrays() + return float(self._solution_col_dual[variable._solver_index]) + + def _constraint_primal(self, constraint): + if not self._has_solution: + return None + self._ensure_solution_arrays() + return float(self._solution_row_value[constraint._solver_index]) + + def _constraint_dual(self, constraint): + if not self._has_solution: + return None + self._ensure_solution_arrays() + return float(self._solution_row_dual[constraint._solver_index]) + + def _get_objective_value(self): + if not self._has_solution: + return None + if self._objective_value is None: + info = self.problem.getInfo() + self._objective_value = float(info.objective_function_value) + return self._objective_value + + +__all__ = ["Variable", "Constraint", "Objective", "Configuration", "Model", "_get_quadratic_terms"] diff --git a/src/optlang/tests/test_highs_interface.py b/src/optlang/tests/test_highs_interface.py new file mode 100644 index 0000000..fa44a9b --- /dev/null +++ b/src/optlang/tests/test_highs_interface.py @@ -0,0 +1,536 @@ +""" +Test suite for the HiGHS solver interface (highs_interface.py), built on +top of optlang's generic abstract test suite (abstract_test_cases.py). + +Scope +----- +highs_interface currently only supports continuous LP/QP problems (see the +module docstring in highs_interface.py) - there is no MILP support. The +abstract suite includes a number of tests that exercise integer/binary +variables and indicator constraints; those are skipped here rather than +implemented, per the LP/QP-only scope of this interface. Everything else +from the abstract suite is implemented. + +`test_clone_model_with_lp` is skipped too: it exercises Model.clone's +use_lp=True path (round-tripping through to_lp()/from_lp()), which this +interface does not implement. + +A couple of the abstract tests are deliberately narrowed for this +interface's design (documented inline where relevant): + - test_tolerance_parameters iterates over dir(tolerances) in the generic + version, which also picks up non-numeric members (methods, dunders); + here it's pinned to the three real tolerance parameters. + - test_config_gets_copied_too checks verbosity/timeout/presolve get + copied by Model.clone, but stops short of tolerances: Configuration's + tolerances property only accepts an actual Tolerances instance, so the + dict Configuration.clone builds for it is silently dropped by the + constructor - a pre-existing gap in clone's tolerances handling that is + out of scope for this change. + +Note on imports: this assumes highs_interface.py sits next to this file (or +on the path) as a standalone module, matching the other project files it +was reviewed alongside. If/when it's moved into the optlang package proper +(e.g. optlang/highs_interface.py), swap the `import highs_interface` line +for `from optlang import highs_interface`. +""" +import copy +import json +import pickle +import unittest + +import highspy + +from optlang import interface, highs_interface + +from optlang.tests import abstract_test_cases + +class HighsVariableTestCase(abstract_test_cases.AbstractVariableTestCase): + interface = highs_interface + + def test_get_primal(self): + self.assertIsNone(self.var.primal) + self.var.lb = 3 + self.var.ub = 10 + self.model.add(self.var) + self.model.objective = self.interface.Objective(self.var, direction="max") + status = self.model.optimize() + self.assertEqual(status, interface.OPTIMAL) + self.assertAlmostEqual(self.var.primal, 10) + + def test_changing_variable_names_is_reflected_in_the_solver(self): + self.model.add(self.var) + self.model.update() + self.var.name = "test_2" + self.assertEqual( + self.model.problem.variableName(self.var._solver_index), "test_2" + ) + self.assertIn("test_2", self.model.variables) + self.assertNotIn("test", self.model.variables) + + def test_set_wrong_type_raises(self): + # Trimmed from the abstract version: this interface is LP/QP only + # (see highs_interface.Variable.__init__), so there's no valid + # "integer" type to switch to for the final two lines of the + # original test. See test_integer_and_binary_types_are_rejected + # below for the MIP-specific replacement. + self.assertRaises(ValueError, self.interface.Variable, name="test", type="mayo") + self.assertRaises(Exception, setattr, self.var, 'type', 'ketchup') + self.model.add(self.var) + self.assertRaises(ValueError, setattr, self.var, "type", "mustard") + + @unittest.skip("HiGHS interface only supports continuous variables (LP/QP, no MIP).") + def test_change_type(self): + pass + + def test_integer_and_binary_types_are_rejected(self): + """Replacement for test_change_type: rather than silently accepting + a variable type the solver can never honor, construction should + fail outright for non-continuous types.""" + self.assertRaises(ValueError, self.interface.Variable, "int_var", type="integer") + self.assertRaises(ValueError, self.interface.Variable, "bin_var", type="binary") + + +class HighsConstraintTestCase(abstract_test_cases.AbstractConstraintTestCase): + interface = highs_interface + + def test_get_primal(self): + self.assertIsNone(self.model.constraints[0].primal) + self.model.optimize() + self.assertEqual(self.model.status, interface.OPTIMAL) + self.assertTrue(isinstance(self.model.constraints[0].primal, float)) + + @unittest.skip( + "Indicator constraints require a binary indicator variable, which " + "this LP/QP-only interface does not support." + ) + def test_indicator_constraint_support(self): + pass + + +class HighsObjectiveTestCase(abstract_test_cases.AbstractObjectiveTestCase): + interface = highs_interface + + def setUp(self): + with open(abstract_test_cases.TESTMODELPATH) as infile: + self.model = self.interface.Model.from_json(json.load(infile)) + self.obj = self.model.objective + + def test_change_direction(self): + self.obj.direction = "min" + self.assertEqual(self.obj.direction, "min") + _, sense = self.model.problem.getObjective() + self.assertEqual(sense, highspy.ObjSense.kMinimize) + + self.obj.direction = "max" + self.assertEqual(self.obj.direction, "max") + _, sense = self.model.problem.getObjective() + self.assertEqual(sense, highspy.ObjSense.kMaximize) + + +class HighsModelTestCase(abstract_test_cases.AbstractModelTestCase): + interface = highs_interface + + # --- abstract methods --------------------------------------------- + + def test_pickle_ability(self): + self.model.optimize() + value = self.model.objective.value + pickle_string = pickle.dumps(self.model) + from_pickle = pickle.loads(pickle_string) + from_pickle.optimize() + self.assertAlmostEqual(from_pickle.objective.value, value) + self.assertEqual( + [(var.lb, var.ub, var.name, var.type) for var in from_pickle.variables.values()], + [(var.lb, var.ub, var.name, var.type) for var in self.model.variables.values()], + ) + self.assertEqual( + [(constr.lb, constr.ub, constr.name) for constr in from_pickle.constraints], + [(constr.lb, constr.ub, constr.name) for constr in self.model.constraints], + ) + + def test_config_gets_copied_too(self): + self.model.configuration.verbosity = 3 + self.model.configuration.timeout = 5 + self.model.configuration.presolve = False + + model_copy = self.interface.Model.clone(self.model) + + self.assertEqual(model_copy.configuration.verbosity, 3) + self.assertEqual(model_copy.configuration.timeout, 5) + self.assertEqual(model_copy.configuration.presolve, False) + + def test_init_from_existing_problem(self): + inner_problem = self.model.problem + self.assertEqual(inner_problem.getNumCol(), len(self.model.variables)) + + new_model = self.interface.Model(problem=inner_problem) + self.assertEqual(len(new_model.variables), len(self.model.variables)) + self.assertEqual(len(new_model.constraints), len(self.model.constraints)) + self.assertEqual( + sorted(var.name for var in new_model.variables), + sorted(var.name for var in self.model.variables), + ) + self.assertEqual( + sorted(c.name for c in new_model.constraints), + sorted(c.name for c in self.model.constraints), + ) + + status = new_model.optimize() + self.assertEqual(status, interface.OPTIMAL) + + def test_add_non_cplex_conform_variable(self): + var = self.interface.Variable('12x!!@#5_3', lb=-666, ub=666) + self.model.add(var) + self.model.update() + self.assertEqual(var.name, self.model.problem.variableName(var._solver_index)) + repickled = pickle.loads(pickle.dumps(self.model)) + var_from_pickle = repickled.variables['12x!!@#5_3'] + self.assertEqual( + var_from_pickle.name, + repickled.problem.variableName(var_from_pickle._solver_index), + ) + + def test_change_of_constraint_is_reflected_in_low_level_solver(self): + x = self.interface.Variable('x', lb=0, ub=10) + y = self.interface.Variable('y', lb=0, ub=10) + z = self.interface.Variable('z', lb=0, ub=10) + constraint = self.interface.Constraint(x + y, lb=0, ub=100, name='constr_test') + self.model.add([z, constraint]) + self.model.update() + + _, idx, _ = self.model.problem.getRowEntries(constraint._solver_index) + self.assertEqual(len(idx), 2) + + constraint.set_linear_coefficients({z: 77.}) + _, idx, val = self.model.problem.getRowEntries(constraint._solver_index) + coeffs = dict(zip( + (self.model.problem.variableName(int(i)) for i in idx), val + )) + self.assertAlmostEqual(coeffs['z'], 77.) + + def test_constraint_set_problem_to_None_caches_the_latest_expression_from_solver_instance(self): + x = self.interface.Variable('x', lb=-83.3, ub=1324422.) + y = self.interface.Variable('y', lb=-181133.3, ub=12000.) + constraint = self.interface.Constraint(0.3 * x + 0.4 * y, lb=-100, name='test') + self.model.add(constraint) + z = self.interface.Variable('z', lb=2, ub=5) + constraint += 77. * z + self.model.remove(constraint) + self.assertEqual( + (constraint.expression - (0.4 * y + 0.3 * x + 77.0 * z)).expand() - 0, 0 + ) + self.assertEqual(constraint.lb, -100) + self.assertEqual(constraint.ub, None) + + def test_change_of_objective_is_reflected_in_low_level_solver(self): + x = self.interface.Variable('x', lb=0, ub=10) + y = self.interface.Variable('y', lb=0, ub=10) + self.model.add([x, y]) + self.model.objective = self.interface.Objective(x + 2 * y, direction='max') + + def _cost_dict(): + obj_expr, _ = self.model.problem.getObjective() + return dict(zip( + (self.model.problem.variableName(int(i)) for i in obj_expr.idxs), + obj_expr.vals, + )) + + coeffs = _cost_dict() + self.assertAlmostEqual(coeffs.get('x', 0), 1.0) + self.assertAlmostEqual(coeffs.get('y', 0), 2.0) + + self.model.objective.set_linear_coefficients({x: 5.0}) + coeffs = _cost_dict() + self.assertAlmostEqual(coeffs.get('x', 0), 5.0) + self.assertAlmostEqual(coeffs.get('y', 0), 2.0) + + def test_change_variable_bounds(self): + inner_problem = self.model.problem + bounds = [(-1000. - i, 1000. + i) for i in range(len(self.model.variables))] + for var, (lb, ub) in zip(self.model.variables, bounds): + # set_bounds updates both at once, avoiding a transient + # lb > (old) ub state that the individual setters would reject. + var.set_bounds(lb, ub) + + lp = inner_problem.getLp() + for i, (lb, ub) in enumerate(bounds): + self.assertAlmostEqual(lp.col_lower_[i], lb) + self.assertAlmostEqual(lp.col_upper_[i], ub) + + def test_change_constraint_bounds(self): + constraint = self.model.constraints[0] + constraint.lb = -1000. + constraint.ub = 1000. + + _, row_lb, row_ub, _ = self.model.problem.getRow(constraint._solver_index) + self.assertAlmostEqual(row_lb, -1000. - constraint._constant) + self.assertAlmostEqual(row_ub, 1000. - constraint._constant) + self.assertAlmostEqual(constraint.lb, -1000.) + self.assertAlmostEqual(constraint.ub, 1000.) + + def test_iadd_objective(self): + x = self.interface.Variable('x', lb=0, ub=10) + y = self.interface.Variable('y', lb=0, ub=10) + self.model.add([x, y]) + + obj = self.interface.Objective(x) + obj += 2. * y + self.model.objective = obj + + obj_expr, _ = self.model.problem.getObjective() + coeffs = dict(zip( + (self.model.problem.variableName(int(i)) for i in obj_expr.idxs), + obj_expr.vals, + )) + self.assertAlmostEqual(coeffs.get('x', 0), 1.0) + self.assertAlmostEqual(coeffs.get('y', 0), 2.0) + + def test_imul_objective(self): + x = self.interface.Variable('x', lb=0, ub=10) + y = self.interface.Variable('y', lb=0, ub=10) + self.model.add([x, y]) + + obj = self.interface.Objective(x + y) + obj *= 3. + self.model.objective = obj + + obj_expr, _ = self.model.problem.getObjective() + coeffs = dict(zip( + (self.model.problem.variableName(int(i)) for i in obj_expr.idxs), + obj_expr.vals, + )) + self.assertAlmostEqual(coeffs.get('x', 0), 3.0) + self.assertAlmostEqual(coeffs.get('y', 0), 3.0) + + def test_set_copied_objective(self): + x = self.interface.Variable('x', lb=0, ub=10) + self.model.add(x) + self.model.objective = self.interface.Objective(x, direction='max') + + obj_copy = copy.copy(self.model.objective) + self.model.objective = obj_copy + + status = self.model.optimize() + self.assertEqual(status, interface.OPTIMAL) + self.assertAlmostEqual(x.primal, 10) + + def test_timeout(self): + self.model.configuration.timeout = 10 + self.assertEqual(self.model.configuration.timeout, 10) + status = self.model.optimize() + self.assertIn(status, (interface.OPTIMAL, interface.TIME_LIMIT)) + + def test_set_linear_coefficients_objective(self): + x = self.interface.Variable('x', lb=0, ub=10) + y = self.interface.Variable('y', lb=0, ub=10) + self.model.add([x, y]) + self.model.objective = self.interface.Objective(x, direction='max') + + self.model.objective.set_linear_coefficients({y: 1.}) + + obj_expr, _ = self.model.problem.getObjective() + coeffs = dict(zip( + (self.model.problem.variableName(int(i)) for i in obj_expr.idxs), + obj_expr.vals, + )) + self.assertAlmostEqual(coeffs.get('x', 0), 1.0) + self.assertAlmostEqual(coeffs.get('y', 0), 1.0) + + def test_set_linear_coefficients_constraint(self): + x = self.interface.Variable('x', lb=0, ub=10) + y = self.interface.Variable('y', lb=0, ub=10) + c = self.interface.Constraint(x, lb=0, ub=100, name='c') + self.model.add([y, c]) + self.model.update() + + c.set_linear_coefficients({y: 3.}) + + _, idx, val = self.model.problem.getRowEntries(c._solver_index) + coeffs = dict(zip( + (self.model.problem.variableName(int(i)) for i in idx), val + )) + self.assertAlmostEqual(coeffs.get('y', 0), 3.0) + + # --- concrete (non-abstract) tests from the base suite that needed + # adapting for this interface: binary/integer variables replaced with + # continuous ones (this interface is LP/QP only), and sloppy=False + # added where a test relies on validation/canonicalization this + # interface skips by default. ----------------------------------- + + def test_add_constraints(self): + x = self.interface.Variable('x', lb=0, ub=1) + y = self.interface.Variable('y', lb=-181133.3, ub=12000.) + z = self.interface.Variable('z', lb=0., ub=3) + constr1 = self.interface.Constraint(0.3 * x + 0.4 * y + 66. * z, lb=-100, ub=0., name='test') + constr2 = self.interface.Constraint(2.333 * x + y + 3.333, ub=100.33, name='test2') + constr3 = self.interface.Constraint(2.333 * x + y + z, ub=100.33, lb=-300) + constr4 = self.interface.Constraint(77 * x, lb=10, name='Mul_constraint') + constr5 = self.interface.Constraint(x, ub=-10, name='Only_var_constraint') + constr6 = self.interface.Constraint(3, ub=88., name='Number_constraint') + self.model.add(constr1) + self.model.add(constr2) + self.model.add(constr3, sloppy=True) + self.model.add([constr4, constr5, constr6]) + self.assertIn(constr1.name, self.model.constraints) + self.assertIn(constr2.name, self.model.constraints) + self.assertIn(constr3.name, self.model.constraints) + self.assertIn(constr4.name, self.model.constraints) + self.assertIn(constr5.name, self.model.constraints) + self.assertIn(constr6.name, self.model.constraints) + + def test_remove_constraints(self): + x = self.interface.Variable('x', lb=0, ub=1) + y = self.interface.Variable('y', lb=-181133.3, ub=12000.) + z = self.interface.Variable('z', lb=4, ub=4) + constr1 = self.interface.Constraint(0.3 * x + 0.4 * y + 66. * z, lb=-100, ub=0., name='test') + self.assertEqual(constr1.problem, None) + self.model.add(constr1) + self.model.update() + self.assertEqual(constr1.problem, self.model) + self.assertIn(constr1, self.model.constraints) + self.model.remove(constr1.name) + self.model.update() + self.assertEqual(constr1.problem, None) + self.assertNotIn(constr1, self.model.constraints) + + def test_add_nonlinear_constraint_raises(self): + x = self.interface.Variable('x', lb=0, ub=1) + y = self.interface.Variable('y', lb=-181133.3, ub=12000.) + z = self.interface.Variable('z', lb=3, ub=3) + with self.assertRaises(ValueError): + constraint = self.interface.Constraint( + 0.3 * x + 0.4 * y ** x + 66. * z, lb=-100, ub=0., name='test' + ) + self.model.add(constraint) + + # --- MIP-only tests from the abstract suite: skipped, see module + # docstring. ------------------------------------------------------ + + @unittest.skip("MIP (integer variables) is not supported by this interface.") + def test_add_integer_var(self): + pass + + @unittest.skip("MIP (integer variables) is not supported by this interface.") + def test_change_variable_type(self): + pass + + @unittest.skip("MIP (binary variables) is not supported by this interface.") + def test_binary_variables(self): + pass + + @unittest.skip("MIP (integer variables) is not supported by this interface.") + def test_integer_variable_dual(self): + pass + + @unittest.skip("MIP (integer variables) is not supported by this interface.") + def test_integer_constraint_dual(self): + pass + + @unittest.skip("MIP (integer variables) is not supported by this interface.") + def test_integer_batch_duals(self): + pass + + @unittest.skip("MIP (integer variables) is not supported by this interface.") + def test_implicitly_convert_milp_to_lp(self): + pass + + @unittest.skip("MIP (integer variables) is not supported by this interface.") + def test_optimize_milp(self): + pass + + @unittest.skip("MIP (integer variables) is not supported by this interface.") + def test_is_integer(self): + pass + + @unittest.skip("Cloning via LP export/import (to_lp/from_lp) is not implemented for this interface.") + def test_clone_model_with_lp(self): + pass + + +class HighsConfigurationTestCase(abstract_test_cases.AbstractConfigurationTestCase): + interface = highs_interface + + def test_tolerance_parameters(self): + # Narrowed from the abstract version, which iterates over + # dir(tolerances) - that also picks up non-numeric members (e.g. + # to_dict, dunders), which "setattr(param, 2 * val)" can't handle. + # Restricting to the actual tunable parameters keeps the same + # intent (round-trip every tolerance through get/set) without that. + model = self.interface.Model() + params = ["feasibility", "optimality", "ipm_optimality"] + for param in params: + val = getattr(model.configuration.tolerances, param) + setattr(model.configuration.tolerances, param, 2 * val) + self.assertEqual(getattr(model.configuration.tolerances, param), 2 * val) + + +class HighsQuadraticProgrammingTestCase(abstract_test_cases.AbstractQuadraticProgrammingTestCase): + interface = highs_interface + + def setUp(self): + self.model = self.interface.Model() + + def test_convex_obj(self): + x = self.interface.Variable('x', lb=-10, ub=10) + y = self.interface.Variable('y', lb=-10, ub=10) + self.model.add([x, y]) + self.model.objective = self.interface.Objective(x ** 2 + y ** 2, direction='min') + + status = self.model.optimize() + self.assertEqual(status, interface.OPTIMAL) + self.assertAlmostEqual(x.primal, 0.0, places=5) + self.assertAlmostEqual(y.primal, 0.0, places=5) + self.assertAlmostEqual(self.model.objective.value, 0.0, places=5) + + def test_non_convex_obj(self): + x = self.interface.Variable('x', lb=-10, ub=10) + y = self.interface.Variable('y', lb=-10, ub=10) + self.model.add([x, y]) + # x^2 - y^2 is indefinite (a saddle at the origin), so this + # objective is not convex. + self.model.objective = self.interface.Objective(x ** 2 - y ** 2, direction='min') + + # HiGHS's QP solver is built for convex objectives. For a + # non-convex one, the exact termination status is solver-version + # dependent (it may converge to a stationary point and still + # report "optimal", report a different status, or raise). We only + # check that solving a non-convex QP doesn't crash the interface + # itself and, if it does return, that the status is one optlang + # recognizes. + try: + status = self.model.optimize() + except Exception: + return + self.assertIn(status, list(interface.statuses.keys())) + + def test_qp_convex(self): + x = self.interface.Variable('x', lb=-10, ub=10) + y = self.interface.Variable('y', lb=-10, ub=10) + self.model.add([x, y]) + self.model.add(self.interface.Constraint(x + y, lb=1, name='c')) + self.model.objective = self.interface.Objective(x ** 2 + y ** 2, direction='min') + + status = self.model.optimize() + self.assertEqual(status, interface.OPTIMAL) + self.assertAlmostEqual(x.primal, 0.5, places=4) + self.assertAlmostEqual(y.primal, 0.5, places=4) + self.assertAlmostEqual(self.model.objective.value, 0.5, places=4) + + def test_qp_non_convex(self): + x = self.interface.Variable('x', lb=-5, ub=5) + y = self.interface.Variable('y', lb=-5, ub=5) + self.model.add([x, y]) + self.model.add(self.interface.Constraint(x + y, lb=-5, ub=5, name='c')) + self.model.objective = self.interface.Objective(x * y, direction='min') + + # See test_non_convex_obj - global optimality isn't guaranteed here, + # so this only checks the solver handles it gracefully. + try: + status = self.model.optimize() + except Exception: + return + self.assertIn(status, list(interface.statuses.keys())) + + +if __name__ == '__main__': + unittest.main() diff --git a/tox.ini b/tox.ini index f78ee40..aa4d152 100644 --- a/tox.ini +++ b/tox.ini @@ -33,7 +33,7 @@ deps = numpy >=1.13 osqp >=1.1.0 mip >=1.9.1; python_version < "3.14" - highspy >=1.5.3 + highspy >=1.14.0 cplex gurobipy symengine: symengine