From ce0328ce365ed97455e953830993907842834547 Mon Sep 17 00:00:00 2001 From: Christian Diener Date: Tue, 23 May 2023 10:48:42 -0700 Subject: [PATCH 01/10] first commit --- src/optlang/osqp_interface.py | 23 +++++++++++++++++++++-- 1 file changed, 21 insertions(+), 2 deletions(-) diff --git a/src/optlang/osqp_interface.py b/src/optlang/osqp_interface.py index 48478bd..1eaf975 100644 --- a/src/optlang/osqp_interface.py +++ b/src/optlang/osqp_interface.py @@ -28,7 +28,8 @@ import numpy as np from numpy import array, concatenate, Infinity, zeros, isnan -from optlang import interface, symbolics, available_solvers +from optlang import symbolics, available_solvers +import optlang.matrix_interface as mi from optlang.util import inheritdocstring from optlang.expression_parsing import parse_optimization_expression from optlang.exceptions import SolverError @@ -77,7 +78,7 @@ _TYPES = ("continuous",) -class OSQPProblem(object): +class OSQPProblem(mi.MatrixProblem): """A concise representation of an OSQP problem. OSQP assumes that the problem will be pretty much immutable. This is @@ -120,6 +121,24 @@ def __init__(self): } self.__solution = None + def map_settings(self): + """Map internal settings to OSQP settings.""" + settings = { + "linsys_solver": "qdldl", + "max_iter": 100000, + "eps_abs": self.settings["optimality_tolerance"], + "eps_rel": self.settings["optimality_tolerance"], + "eps_prim_inf": self.settings["primal_inf_tolerance"], + "eps_dual_inf": self.settings["dual_inf_tolerance"], + "polish": True, + "verbose": False, + "scaling": 0, + "time_limit": 0, + "adaptive_rho": True, + "rho": 1.0, + "alpha": 1.6, + } + def build(self): """Build the problem instance.""" d = self.direction From d22957226d0514d6ccc7e87383fb099330a50efd Mon Sep 17 00:00:00 2001 From: Christian Diener Date: Tue, 23 May 2023 10:48:53 -0700 Subject: [PATCH 02/10] first commit --- src/optlang/matrix_interface.py | 896 ++++++++++++++++++++++++++++++++ 1 file changed, 896 insertions(+) create mode 100644 src/optlang/matrix_interface.py diff --git a/src/optlang/matrix_interface.py b/src/optlang/matrix_interface.py new file mode 100644 index 0000000..9429112 --- /dev/null +++ b/src/optlang/matrix_interface.py @@ -0,0 +1,896 @@ +# Copyright 2013 Novo Nordisk Foundation Center for Biosustainability, +# Technical University of Denmark. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Generic solver interface for (sparse) matrix solvers. + +This interface is only supposed to be a basis to add matrix solvers to optlang. It +creates a thin layer between all optlang operations and the standard formulation of +LP, QP, and MILP problems. +""" + +import logging +import six +import pickle +import numpy as np +from numpy import array, concatenate, Infinity, zeros + +from optlang import interface, symbolics, available_solvers +from optlang.util import inheritdocstring +from optlang.expression_parsing import parse_optimization_expression +from optlang.exceptions import SolverError + +from scipy.sparse import csc_matrix + +from optlang.symbolics import add, mul + +log = logging.getLogger(__name__) + + +_STATUS_MAP = { + "interrupted by user": interface.ABORTED, + "run time limit reached": interface.TIME_LIMIT, + "feasible": interface.FEASIBLE, + "primal infeasible": interface.INFEASIBLE, + "dual infeasible": interface.INFEASIBLE, + "primal infeasible inaccurate": interface.INFEASIBLE, + "dual infeasible inaccurate": interface.INFEASIBLE, + "solved inaccurate": interface.NUMERIC, + "solved": interface.OPTIMAL, + "maximum iterations reached": interface.ITERATION_LIMIT, + "unsolved": interface.SPECIAL, + "problem non convex": interface.SPECIAL, + "non-existing-status": "Here for testing that missing statuses are handled.", +} + +_LP_METHODS = ("auto", ) + +_QP_METHODS = ("auto", ) + +_TYPES = ("continuous", "binary", "integer") + + +class MatrixProblem(object): + """A concise representation of an LP/QP/MILP problem. + + This class assumes that the problem will be pretty much immutable. This is + a small intermediate layer based on dictionaries that is fast to modify + but can also be converted to a matrix formulation quickly. This is the only + thing you will need to modify. + """ + + def __init__(self): + self.variables = set() + self.integer_vars = set() + self.constraints = set() + self.constraint_coefs = dict() + self.constraint_lbs = dict() + self.constraint_ubs = dict() + self.variable_lbs = dict() + self.variable_ubs = dict() + self.obj_linear_coefs = dict() + self.obj_quadratic_coefs = dict() + self.primals = {} + self.cprimals = {} + self.duals = {} + self.vduals = {} + self.obj_value = None + self.direction = -1 + self.info = None + self.status = None + self.__solution = None + self.default_settings() + + def default_settings(self): + """Set the default settings of the solver.""" + self.settings = { + "presolve": True, + "lp_method": "auto", + "qp_method": "auto", + "time_limit": 0, + "primal_inf_tolerance": 1e-6, + "dual_inf_tolerance": 1e-6, + "optimality_tolerance": 1e-6, + "verbose": 0 + } + + def build(self): + """Build the problem instance.""" + d = self.direction + vmap = dict(zip(self.variables, range(len(self.variables)))) + nv = len(self.variables) + cmap = dict(zip(self.constraints, range(len(self.constraints)))) + nc = len(self.constraints) + if len(self.obj_quadratic_coefs) > 0: + P = array(raise NotImplementedError("This needs to be overwritten by the child class.") + [ + [vmap[vn[0]], vmap[vn[1]], coef * d * 2.0] + for vn, coef in six.iteritems(self.obj_quadratic_coefs) + ] + ) + P = csc_matrix( + (P[:, 2], (P[:, 0].astype("int64"), P[:, 1].astype("int64"))), + shape=(nv, nv), + ) + else: + P = None + q = zeros(nv) + q[[vmap[vn] for vn in self.obj_linear_coefs]] = list( + self.obj_linear_coefs.values() + ) + q = q * d + Av = array([[vmap[k] + nc, vmap[k], 1.0] for k in self.variables]) + vbounds = array( + [[self.variable_lbs[vn], self.variable_ubs[vn]] for vn in self.variables] + ) + if len(self.constraint_coefs) > 0: + A = array( + [ + [cmap[vn[0]], vmap[vn[1]], coef] + for vn, coef in six.iteritems(self.constraint_coefs) + ] + ) + bounds = array( + [ + [self.constraint_lbs[cn], self.constraint_ubs[cn]] + for cn in self.constraints + ] + ) + A = concatenate((A, Av)) + bounds = concatenate((bounds, vbounds)) + else: + A = Av + bounds = vbounds + if A.shape[0] == 0: + A = None + else: + A = csc_matrix( + (A[:, 2], (A[:, 0].astype("int64"), A[:, 1].astype("int64"))), + shape=(nc + nv, nv), + ) + return P, q, A, bounds + + def solve(self): + """Solve the problem.""" + raise NotImplementedError("This needs to be overwritten by the child class.") + + def reset(self, everything=False): + """Reset the public solver solution.""" + self.info = None + self.primals = {} + self.cprimals = {} + self.duals = {} + self.vduals = {} + if everything: + self.__solution = None + + def still_valid(self, A, bounds): + """Check if previous solutions is still feasible.""" + if len(self.__solution["x"]) != len(self.variables) or len( + self.__solution["y"] + ) != len(self.constraints): + return False + c = A.dot(self.__solution["x"]) + ea = self.settings["eps_abs"] + er = self.settings["eps_rel"] + valid = np.all( + (c + er * np.abs(c) + ea >= bounds[:, 0]) + & (c - er * np.abs(c) - ea <= bounds[:, 1]) + ) + return valid + + def clean(self): + """Remove unused variables and constraints.""" + self.reset() + self.variable_lbs = { + k: v for k, v in six.iteritems(self.variable_lbs) if k in self.variables + } + self.variable_ubs = { + k: v for k, v in six.iteritems(self.variable_ubs) if k in self.variables + } + self.constraint_coefs = { + k: v + for k, v in six.iteritems(self.constraint_coefs) + if k[0] in self.constraints and k[1] in self.variables + } + self.obj_linear_coefs = { + k: v for k, v in six.iteritems(self.obj_linear_coefs) if k in self.variables + } + self.obj_quadratic_coefs = { + k: v + for k, v in six.iteritems(self.obj_quadratic_coefs) + if k[0] in self.variables and k[1] in self.variables + } + self.constraint_lbs = { + k: v for k, v in six.iteritems(self.constraint_lbs) if k in self.constraints + } + self.constraint_ubs = { + k: v for k, v in six.iteritems(self.constraint_ubs) if k in self.constraints + } + + def rename_constraint(self, old, new): + """Rename a constraint.""" + self.reset() + self.constraints.remove(old) + self.constraints.add(new) + name_map = {k: k for k in self.constraints} + name_map[old] = new + self.constraint_coefs = { + (name_map[k[0]], k[1]): v for k, v in six.iteritems(self.constraint_coefs) + } + self.constraint_lbs[new] = self.constraint_lbs.pop(old) + self.constraint_ubs[new] = self.constraint_ubs.pop(old) + + def rename_variable(self, old, new): + """Rename a variable.""" + self.reset() + self.variables.remove(old) + self.variables.add(new) + name_map = {k: k for k in self.variables} + name_map[old] = new + self.constraint_coefs = { + (k[0], name_map[k[1]]): v for k, v in six.iteritems(self.constraint_coefs) + } + self.obj_quadratic_coefs = { + (name_map[k[0]], name_map[k[1]]): v + for k, v in six.iteritems(self.obj_quadratic_coefs) + } + self.obj_linear_coefs = { + name_map[k]: v for k, v in six.iteritems(self.obj_linear_coefs) + } + self.variable_lbs[new] = self.variable_lbs.pop(old) + self.variable_ubs[new] = self.variable_ubs.pop(old) + + +@six.add_metaclass(inheritdocstring) +class Variable(interface.Variable): + def __init__(self, name, *args, **kwargs): + super(Variable, self).__init__(name, **kwargs) + + @interface.Variable.type.setter + def type(self, value): + if self.problem is not None: + if value not in _TYPES: + raise ValueError( + "This interface cannot handle variables of type '%s'. " % value + + "The following variable types are available: " + + ", ".join(_TYPES) + ) + super(Variable, Variable).type.fset(self, value) + + def _get_primal(self): + return self.problem.problem.primals.get(self.name, None) + + @property + def dual(self): + if self.problem is not None: + return self.problem.problem.vduals.get(self.name, None) + return None + + @interface.Variable.name.setter + def name(self, value): + old_name = getattr(self, "name", None) + super(Variable, Variable).name.fset(self, value) + if getattr(self, "problem", None) is not None: + if old_name != value: + self.problem.problem.rename_variable(old_name, value) + + +@six.add_metaclass(inheritdocstring) +class Constraint(interface.Constraint): + _INDICATOR_CONSTRAINT_SUPPORT = False + + def __init__(self, expression, sloppy=False, *args, **kwargs): + super(Constraint, self).__init__(expression, *args, sloppy=sloppy, **kwargs) + + def set_linear_coefficients(self, coefficients): + if self.problem is not None: + self.problem.update() + self.problem.problem.reset() + for var, coef in six.iteritems(coefficients): + self.problem.problem.constraint_coefs[(self.name, var.name)] = float( + coef + ) + else: + raise Exception( + "Can't change coefficients if constraint is not associated " + "with a model." + ) + + def get_linear_coefficients(self, variables): + if self.problem is not None: + self.problem.update() + coefs = { + v: self.problem.problem.constraint_coefs.get((self.name, v.name), 0.0) + for v in variables + } + return coefs + else: + raise Exception( + "Can't get coefficients from solver if constraint is not " "in a model" + ) + + def _get_expression(self): + if self.problem is not None: + variables = self.problem._variables + all_coefs = self.problem.problem.constraint_coefs + coefs = [(v, all_coefs.get((self.name, v.name), 0.0)) for v in variables] + expression = add([mul((symbolics.Real(co), v)) for (v, co) in coefs]) + self._expression = expression + return self._expression + + @property + def problem(self): + return self._problem + + @problem.setter + def problem(self, value): + if value is None: + # Update expression from solver instance one last time + self._get_expression() + self._problem = value + + @property + def primal(self): + if self.problem is None: + return None + return self.problem.problem.cprimals.get(self.name, None) + + @property + def dual(self): + if self.problem is None: + return None + d = self.problem.problem.duals.get(self.name, None) + if d is not None: + d = -d + return d + + @interface.Constraint.name.setter + def name(self, value): + old_name = self.name + super(Constraint, Constraint).name.fset(self, value) + if getattr(self, "problem", None) is not None: + if old_name != value: + self.problem.problem.rename_constraint(old_name, value) + + @interface.Constraint.lb.setter + def lb(self, value): + self._check_valid_lower_bound(value) + if getattr(self, "problem", None) is not None: + lb = -Infinity if value is None else float(value) + self.problem.problem.constraint_lbs[self.name] = lb + self._lb = value + + @interface.Constraint.ub.setter + def ub(self, value): + self._check_valid_upper_bound(value) + if getattr(self, "problem", None) is not None: + ub = Infinity if value is None else float(value) + self.problem.problem.constraint_ubs[self.name] = ub + self._ub = value + + def __iadd__(self, other): + # if self.problem is not None: + # self.problem._add_to_constraint(self.index, other) + if self.problem is not None: + problem_reference = self.problem + self.problem._remove_constraint(self) + super(Constraint, self).__iadd__(other) + problem_reference._add_constraint(self, sloppy=False) + else: + super(Constraint, self).__iadd__(other) + return self + + +@six.add_metaclass(inheritdocstring) +class Objective(interface.Objective): + def __init__(self, expression, sloppy=False, **kwargs): + super(Objective, self).__init__(expression, sloppy=sloppy, **kwargs) + self._expression_expired = False + if not (sloppy or self.is_Linear or self.is_Quadratic): + raise ValueError( + "This interface only supports linear and quadratic objectives.") + + @property + def value(self): + if getattr(self, "problem", None) is None: + return None + if self.problem.problem.obj_value is None: + return None + return self.problem.problem.obj_value + self.problem._objective_offset + + @interface.Objective.direction.setter + def direction(self, value): + super(Objective, Objective).direction.__set__(self, value) + if value == "min": + self.problem.problem.direction = 1 + else: + self.problem.problem.direction = -1 + + def _get_expression(self): + if ( + self.problem is not None + and self._expression_expired + and len(self.problem._variables) > 0 + ): + model = self.problem + vars = model._variables + expression = add( + [ + coef * vars[vn] + for vn, coef in six.iteritems(model.problem.obj_linear_coefs) + ] + ) + q_ex = add( + [ + coef * vars[vn[0]] * vars[vn[1]] + for vn, coef in six.iteritems(model.problem.obj_quadratic_coefs) + ] + ) + expression += q_ex + self._expression = expression + self._expression_expired = False + return self._expression + + def set_linear_coefficients(self, coefficients): + if self.problem is not None: + self.problem.update() + for v, coef in six.iteritems(coefficients): + self.problem.problem.obj_linear_coefs[v.name] = float(coef) + self._expression_expired = True + else: + raise Exception( + "Can't change coefficients if objective is not associated " + "with a model." + ) + + def get_linear_coefficients(self, variables): + if self.problem is not None: + self.problem.update() + coefs = { + v: self.problem.problem.obj_linear_coefs.get(v.name, 0.0) + for v in variables + } + return coefs + else: + raise Exception( + "Can't get coefficients from solver if objective " "is not in a model" + ) + + +@six.add_metaclass(inheritdocstring) +class Configuration(interface.MathematicalProgrammingConfiguration): + def __init__( + self, + lp_method="auto", + presolve=False, + verbosity=0, + timeout=None, + qp_method="auto", + *args, + **kwargs + ): + super(Configuration, self).__init__(*args, **kwargs) + self.lp_method = lp_method + self.presolve = presolve + self.verbosity = verbosity + self.timeout = timeout + self.qp_method = qp_method + if "tolerances" in kwargs: + for key, val in six.iteritems(kwargs["tolerances"]): + if key in self._tolerance_functions(): + setattr(self.tolerances, key, val) + + @property + def lp_method(self): + """The algorithm used to solve LP problems.""" + return "default" + + @lp_method.setter + def lp_method(self, lp_method): + if lp_method not in _LP_METHODS: + raise ValueError( + "LP Method %s is not valid (choose one of: %s)" + % (lp_method, ", ".join(_LP_METHODS)) + ) + + @property + def linear_solver(self): + return self._linear_solver + + @linear_solver.setter + def linear_solver(self, solver): + if solver not in ("qdldl", "mkl pardiso"): + raise ValueError( + "%s is not valid (choose either `qdldl` or `mkl pardiso`)" % solver + ) + if getattr(self, "problem", None) is not None: + self.problem.problem.settings["linsys_solver"] = solver + self._linear_solver = solver + + def _set_presolve(self, value): + if getattr(self, "problem", None) is not None: + if value is True: + self.problem.problem.settings["scaling"] = 10 + elif value is False or value == "auto": + self.problem.problem.settings["scaling"] = 0 + else: + raise ValueError( + str(value) + " is not a valid presolve parameter. Must be True, " + "False or 'auto'." + ) + + @property + def presolve(self): + return self.problem.problem.settings["presolve"] + + @presolve.setter + def presolve(self, value): + self._set_presolve(value) + self._presolve = value + + @property + def verbosity(self): + return self._verbosity + + @verbosity.setter + def verbosity(self, value): + if getattr(self, "problem", None) is not None: + self.problem.problem.settings["verbose"] = int(value) + self._verbosity = value + + @property + def timeout(self): + return self.problem.problem.settings["time_limit"] + + @timeout.setter + def timeout(self, value): + if getattr(self, "problem", None) is not None: + if value is None: + self.problem.problem.settings["time_limit"] = 0 + else: + self.problem.problem.settings["time_limit"] = value + + def __getstate__(self): + return { + "presolve": self.presolve, + "timeout": self.timeout, + "verbosity": self.verbosity, + "linear_solver": self.linear_solver, + "tolerances": { + "feasibility": self.tolerances.feasibility, + "optimality": self.tolerances.optimality, + "integrality": self.tolerances.integrality, + }, + } + + def __setstate__(self, state): + for key, val in six.iteritems(state): + if key != "tolerances": + setattr(self, key, val) + for key, val in six.iteritems(state["tolerances"]): + if key in self._tolerance_functions(): + setattr(self.tolerances, key, val) + + @property + def qp_method(self): + """Change the algorithm used to optimize QP problems.""" + return "auto" + + @qp_method.setter + def qp_method(self, value): + if value not in _QP_METHODS: + raise ValueError( + "%s is not a valid qp_method. Choose between %s" + % (value, str(_QP_METHODS)) + ) + self._qp_method = value + + def _get_feasibility(self): + return self.problem.problem.settings["primal_infeasibility"] + + def _set_feasibility(self, value): + self.problem.problem.settings["primal_inf_tolerance"] = value + self.problem.problem.settings["dual_inf_tolerance"] = value + + def _get_integrality(self): + return 1e-6 + + def _set_integrality(self, value): + pass + + def _get_optimality(self): + return self.problem.problem.settings["optimality_tolerance"] + + def _set_optimality(self, value): + self.problem.problem.settings["optimality_tolerance"] = value + + def _tolerance_functions(self): + return { + "feasibility": (self._get_feasibility, self._set_feasibility), + "optimality": (self._get_optimality, self._set_optimality), + "integrality": (self._get_integrality, self._set_integrality), + } + + +@six.add_metaclass(inheritdocstring) +class Model(interface.Model): + ProblemClass = MatrixProblem + + def _initialize_problem(self): + self.problem = self.ProblemClass() + + def _initialize_model_from_problem(self, problem, vc_mapping=None, offset=0): + if not isinstance(problem, self.ProblemClass): + raise TypeError("Provided problem is not a valid OSQP model.") + self.problem = problem + for name in self.problem.variables: + var = Variable( + name, + lb=self.problem.variable_lbs[name], + ub=self.problem.variable_ubs[name], + problem=self, + ) + super(Model, self)._add_variables([var]) + + for name in self.problem.constraints: + # Since constraint expressions are lazily retrieved from the + # solver they don't have to be built here + lhs = symbolics.Integer(0) + constr = Constraint( + lhs, + lb=self.problem.constraint_lbs[name], + ub=self.problem.constraint_ubs[name], + name=name, + problem=self, + ) + + super(Model, self)._add_constraints([constr], sloppy=True) + + if vc_mapping is None: + for constr in self.constraints: + name = constr.name + for variable in constr.variables: + try: + self._variables_to_constraints_mapping[variable.name].add(name) + except KeyError: + self._variables_to_constraints_mapping[variable.name] = set( + [name] + ) + else: + self._variables_to_constraints_mapping = vc_mapping + + linear_expression = add( + [ + coef * self._variables[vn] + for vn, coef in six.iteritems(self.problem.obj_linear_coefs) + ] + ) + quadratic_expression = add( + [ + coef * self._variables[vn[0]] * self._variables[vn[1]] + for vn, coef in six.iteritems(self.problem.obj_quadratic_coefs) + ] + ) + + self._objective_offset = offset + self._objective = Objective( + linear_expression + quadratic_expression + offset, + problem=self, + direction={-1: "max", 1: "min"}[self.problem.direction], + name="osqp_objective", + ) + + @property + def objective(self): + return self._objective + + @objective.setter + def objective(self, value): + value.problem = None + if self._objective is not None: # Reset previous objective + self.problem.obj_linear_coefs = dict() + self.problem.obj_quadratic_coefs = dict() + super(Model, self.__class__).objective.fset(self, value) + self.update() + expression = self._objective._expression + ( + offset, + linear_coefficients, + quadratic_coeffients, + ) = parse_optimization_expression(value, quadratic=True, expression=expression) + self._objective_offset = offset + if linear_coefficients: + self.problem.obj_linear_coefs = { + v.name: float(c) for v, c in six.iteritems(linear_coefficients) + } + + for key, coef in six.iteritems(quadratic_coeffients): + if len(key) == 1: + var = six.next(iter(key)) + self.problem.obj_quadratic_coefs[(var.name, var.name)] = float(coef) + else: + var1, var2 = key + self.problem.obj_quadratic_coefs[(var1.name, var2.name)] = 0.5 * float( + coef + ) + self.problem.obj_quadratic_coefs[(var2.name, var1.name)] = 0.5 * float( + coef + ) + + self._set_objective_direction(value.direction) + value.problem = self + + def _set_objective_direction(self, direction): + self.problem.direction = {"min": 1, "max": -1}[direction] + + def _get_primal_values(self): + if len(self.problem.primals) == 0: + raise SolverError("The problem has not been solved yet!") + primal_values = [self.problem.primals[v.name] for v in self._variables] + return primal_values + + def _get_reduced_costs(self): + if len(self.problem.duals) == 0: + raise SolverError("The problem has not been solved yet!") + reduced_costs = [self.problem.vduals[v.name] for v in self._variables] + return reduced_costs + + def _get_constraint_values(self): + if len(self.problem.primals) == 0: + raise SolverError("The problem has not been solved yet!") + constraint_primals = [self.problem.cprimals[c.name] for c in self._constraints] + return constraint_primals + + def _get_shadow_prices(self): + if len(self.problem.primals) == 0: + raise SolverError("The problem has not been solved yet!") + dual_values = [-self.problem.duals[c.name] for c in self._constraints] + return dual_values + + @property + def is_integer(self): + return False + + def _optimize(self): + self.update() + self.problem.solve() + prior_status = self.problem.status + self._original_status = prior_status + status = _STATUS_MAP[prior_status] + return status + + def _set_variable_bounds_on_problem(self, var_lb, var_ub): + self.problem.reset() + for var, val in var_lb: + lb = -Infinity if val is None else float(val) + self.problem.variable_lbs[var.name] = lb + for var, val in var_ub: + ub = Infinity if val is None else val + self.problem.variable_ubs[var.name] = float(ub) + + def _add_variables(self, variables): + super(Model, self)._add_variables(variables) + self.problem.reset() + for variable in variables: + lb = -Infinity if variable.lb is None else float(variable.lb) + ub = Infinity if variable.ub is None else float(variable.ub) + self.problem.variables.add(variable.name) + self.problem.variable_lbs[variable.name] = lb + self.problem.variable_ubs[variable.name] = ub + variable.problem = self + + def _remove_variables(self, variables): + # Not calling parent method to avoid expensive variable removal + # from sympy expressions + if self.objective is not None: + self.objective._expression = self.objective.expression.xreplace( + {var: 0 for var in variables} + ) + for variable in variables: + del self._variables_to_constraints_mapping[variable.name] + variable.problem = None + del self._variables[variable.name] + self.problem.variables.remove(variable.name) + self.problem.clean() + + def _add_constraints(self, constraints, sloppy=False): + super(Model, self)._add_constraints(constraints, sloppy=sloppy) + self.problem.reset() + for constraint in constraints: + # This needs to be done in order to not trigger constraint._get_expression() + constraint._problem = None + if constraint.is_Linear: + _, coeff_dict, _ = parse_optimization_expression(constraint) + lb = -Infinity if constraint.lb is None else float(constraint.lb) + ub = Infinity if constraint.ub is None else float(constraint.ub) + self.problem.constraints.add(constraint.name) + self.problem.constraint_coefs.update( + { + (constraint.name, v.name): float(co) + for v, co in six.iteritems(coeff_dict) + } + ) + self.problem.constraint_lbs[constraint.name] = lb + self.problem.constraint_ubs[constraint.name] = ub + constraint.problem = self + elif constraint.is_Quadratic: + raise NotImplementedError( + "Quadratic constraints (like %s) are not supported " + "in this interface yet." % constraint + ) + else: + raise ValueError( + "This interface only supports linear or quadratic constraints. " + "%s is neither linear nor quadratic." % constraint + ) + + def _remove_constraints(self, constraints): + super(Model, self)._remove_constraints(constraints) + for constraint in constraints: + self.problem.constraints.remove(constraint.name) + self.problem.clean() + + def _get_variable_indices(self, names): + vmap = dict(zip(self.variables, range(len(self.variables)))) + return [vmap[n] for n in names] + + def __setstate__(self, d): + self.__init__() + problem = pickle.loads(d["problem"]) + # workaround to conserve the order + problem.variables = d["variables"] + problem.constraints = d["constraints"] + self._initialize_model_from_problem( + problem, vc_mapping=d["v_to_c"], offset=d["offset"] + ) + problem.variables = set(problem.variables) + problem.constraints = set(problem.constraints) + self.configuration = Configuration() + self.configuration.problem = self + self.configuration.__setstate__(d["config"]) + + def __getstate__(self): + self.problem.reset() + self.update() + return { + "problem": pickle.dumps(self.problem), + "variables": [v.name for v in self._variables], + "constraints": [c.name for c in self._constraints], + "v_to_c": self._variables_to_constraints_mapping, + "config": self.configuration.__getstate__(), + "offset": getattr(self, "_objective_offset", 0.0), + } + + @classmethod + def from_lp(self, lp_problem_str): + """Read a model from an LP file. + + The solver may not have an integrated LP reader so it will either use + cplex or glpk to read the model. This means that QP problems will + currently require cplex to be read :( + """ + if available_solvers["CPLEX"]: + from optlang import cplex_interface + + mod = cplex_interface.Model.from_lp(lp_problem_str) + mod.configuration.lp_method = "auto" + mod.configuration.qp_method = "auto" + return super(Model, self).clone(mod) + else: + from optlang import glpk_interface + + mod = glpk_interface.Model.from_lp(lp_problem_str) + mod.configuration.lp_method = "auto" + return super(Model, self).clone(mod) From 266ef890236747f328b41bb96d80bcbe1763cd9c Mon Sep 17 00:00:00 2001 From: Christian Diener Date: Wed, 30 Aug 2023 12:13:13 -0700 Subject: [PATCH 03/10] first refactor draft --- src/optlang/matrix_interface.py | 52 +- src/optlang/osqp_interface.py | 857 +---------------------- src/optlang/tests/test_osqp_interface.py | 16 +- 3 files changed, 38 insertions(+), 887 deletions(-) diff --git a/src/optlang/matrix_interface.py b/src/optlang/matrix_interface.py index 9429112..1ab2449 100644 --- a/src/optlang/matrix_interface.py +++ b/src/optlang/matrix_interface.py @@ -1,6 +1,3 @@ -# Copyright 2013 Novo Nordisk Foundation Center for Biosustainability, -# Technical University of Denmark. -# # Licensed under the Apache License, Version 2.0 (the "License"); # you may not use this file except in compliance with the License. # You may obtain a copy of the License at @@ -89,7 +86,7 @@ def __init__(self): self.direction = -1 self.info = None self.status = None - self.__solution = None + self._solution = None self.default_settings() def default_settings(self): @@ -113,7 +110,7 @@ def build(self): cmap = dict(zip(self.constraints, range(len(self.constraints)))) nc = len(self.constraints) if len(self.obj_quadratic_coefs) > 0: - P = array(raise NotImplementedError("This needs to be overwritten by the child class.") + P = array( [ [vmap[vn[0]], vmap[vn[1]], coef * d * 2.0] for vn, coef in six.iteritems(self.obj_quadratic_coefs) @@ -173,22 +170,11 @@ def reset(self, everything=False): self.duals = {} self.vduals = {} if everything: - self.__solution = None + self._solution = None def still_valid(self, A, bounds): """Check if previous solutions is still feasible.""" - if len(self.__solution["x"]) != len(self.variables) or len( - self.__solution["y"] - ) != len(self.constraints): - return False - c = A.dot(self.__solution["x"]) - ea = self.settings["eps_abs"] - er = self.settings["eps_rel"] - valid = np.all( - (c + er * np.abs(c) + ea >= bounds[:, 0]) - & (c - er * np.abs(c) - ea <= bounds[:, 1]) - ) - return valid + raise NotImplementedError("This needs to be overwritten by the child class.") def clean(self): """Remove unused variables and constraints.""" @@ -218,6 +204,7 @@ def clean(self): self.constraint_ubs = { k: v for k, v in six.iteritems(self.constraint_ubs) if k in self.constraints } + self.integer_vars = {v for v in self.integer_vars if v in self.variables} def rename_constraint(self, old, new): """Rename a constraint.""" @@ -495,7 +482,7 @@ def __init__( @property def lp_method(self): """The algorithm used to solve LP problems.""" - return "default" + return "auto" @lp_method.setter def lp_method(self, lp_method): @@ -505,26 +492,12 @@ def lp_method(self, lp_method): % (lp_method, ", ".join(_LP_METHODS)) ) - @property - def linear_solver(self): - return self._linear_solver - - @linear_solver.setter - def linear_solver(self, solver): - if solver not in ("qdldl", "mkl pardiso"): - raise ValueError( - "%s is not valid (choose either `qdldl` or `mkl pardiso`)" % solver - ) - if getattr(self, "problem", None) is not None: - self.problem.problem.settings["linsys_solver"] = solver - self._linear_solver = solver - def _set_presolve(self, value): if getattr(self, "problem", None) is not None: - if value is True: - self.problem.problem.settings["scaling"] = 10 - elif value is False or value == "auto": - self.problem.problem.settings["scaling"] = 0 + if isinstance(value, bool): + self.problem.problem.settings["presolve"] = value + elif value == "auto": + self.problem.problem.settings["presove"] = False else: raise ValueError( str(value) + " is not a valid presolve parameter. Must be True, " @@ -567,7 +540,6 @@ def __getstate__(self): "presolve": self.presolve, "timeout": self.timeout, "verbosity": self.verbosity, - "linear_solver": self.linear_solver, "tolerances": { "feasibility": self.tolerances.feasibility, "optimality": self.tolerances.optimality, @@ -598,7 +570,7 @@ def qp_method(self, value): self._qp_method = value def _get_feasibility(self): - return self.problem.problem.settings["primal_infeasibility"] + return self.problem.problem.settings["primal_inf_tolerance"] def _set_feasibility(self, value): self.problem.problem.settings["primal_inf_tolerance"] = value @@ -857,7 +829,7 @@ def __setstate__(self, d): ) problem.variables = set(problem.variables) problem.constraints = set(problem.constraints) - self.configuration = Configuration() + self.configuration = self.interface.Configuration() self.configuration.problem = self self.configuration.__setstate__(d["config"]) diff --git a/src/optlang/osqp_interface.py b/src/optlang/osqp_interface.py index 1eaf975..4e68033 100644 --- a/src/optlang/osqp_interface.py +++ b/src/optlang/osqp_interface.py @@ -1,6 +1,3 @@ -# Copyright 2013 Novo Nordisk Foundation Center for Biosustainability, -# Technical University of Denmark. -# # Licensed under the Apache License, Version 2.0 (the "License"); # you may not use this file except in compliance with the License. # You may obtain a copy of the License at @@ -24,19 +21,10 @@ """ import logging import six -import pickle import numpy as np -from numpy import array, concatenate, Infinity, zeros, isnan -from optlang import symbolics, available_solvers import optlang.matrix_interface as mi from optlang.util import inheritdocstring -from optlang.expression_parsing import parse_optimization_expression -from optlang.exceptions import SolverError - -from scipy.sparse import csc_matrix - -from optlang.symbolics import add, mul log = logging.getLogger(__name__) @@ -55,29 +43,6 @@ raise ImportError("The osqp_interface requires osqp or cuosqp!") -_STATUS_MAP = { - "interrupted by user": interface.ABORTED, - "run time limit reached": interface.TIME_LIMIT, - "feasible": interface.FEASIBLE, - "primal infeasible": interface.INFEASIBLE, - "dual infeasible": interface.INFEASIBLE, - "primal infeasible inaccurate": interface.INFEASIBLE, - "dual infeasible inaccurate": interface.INFEASIBLE, - "solved inaccurate": interface.NUMERIC, - "solved": interface.OPTIMAL, - "maximum iterations reached": interface.ITERATION_LIMIT, - "unsolved": interface.SPECIAL, - "problem non convex": interface.SPECIAL, - "non-existing-status": "Here for testing that missing statuses are handled.", -} - -_LP_METHODS = ("auto", "primal") - -_QP_METHODS = ("auto", "primal") - -_TYPES = ("continuous",) - - class OSQPProblem(mi.MatrixProblem): """A concise representation of an OSQP problem. @@ -86,41 +51,6 @@ class OSQPProblem(mi.MatrixProblem): but can also be converted to an OSQP problem without too much hassle. """ - def __init__(self): - self.variables = set() - self.constraints = set() - self.constraint_coefs = dict() - self.constraint_lbs = dict() - self.constraint_ubs = dict() - self.variable_lbs = dict() - self.variable_ubs = dict() - self.obj_linear_coefs = dict() - self.obj_quadratic_coefs = dict() - self.primals = {} - self.cprimals = {} - self.duals = {} - self.vduals = {} - self.obj_value = None - self.direction = -1 - self.info = None - self.status = None - self.settings = { - "linsys_solver": "qdldl", - "max_iter": 100000, - "eps_abs": 1e-6, - "eps_rel": 1e-6, - "eps_prim_inf": 1e-6, - "eps_dual_inf": 1e-6, - "polish": True, - "verbose": False, - "scaling": 0, - "time_limit": 0, - "adaptive_rho": True, - "rho": 1.0, - "alpha": 1.6, - } - self.__solution = None - def map_settings(self): """Map internal settings to OSQP settings.""" settings = { @@ -131,83 +61,28 @@ def map_settings(self): "eps_prim_inf": self.settings["primal_inf_tolerance"], "eps_dual_inf": self.settings["dual_inf_tolerance"], "polish": True, - "verbose": False, - "scaling": 0, - "time_limit": 0, + "verbose": int(self.settings["verbose"] > 0), + "scaling": 10 if self.settings["presolve"] else 0, + "time_limit": self.settings["time_limit"], "adaptive_rho": True, "rho": 1.0, "alpha": 1.6, } - - def build(self): - """Build the problem instance.""" - d = self.direction - vmap = dict(zip(self.variables, range(len(self.variables)))) - nv = len(self.variables) - cmap = dict(zip(self.constraints, range(len(self.constraints)))) - nc = len(self.constraints) - if len(self.obj_quadratic_coefs) > 0: - P = array( - [ - [vmap[vn[0]], vmap[vn[1]], coef * d * 2.0] - for vn, coef in six.iteritems(self.obj_quadratic_coefs) - ] - ) - P = csc_matrix( - (P[:, 2], (P[:, 0].astype("int64"), P[:, 1].astype("int64"))), - shape=(nv, nv), - ) - else: - P = None - q = zeros(nv) - q[[vmap[vn] for vn in self.obj_linear_coefs]] = list( - self.obj_linear_coefs.values() - ) - q = q * d - Av = array([[vmap[k] + nc, vmap[k], 1.0] for k in self.variables]) - vbounds = array( - [[self.variable_lbs[vn], self.variable_ubs[vn]] for vn in self.variables] - ) - if len(self.constraint_coefs) > 0: - A = array( - [ - [cmap[vn[0]], vmap[vn[1]], coef] - for vn, coef in six.iteritems(self.constraint_coefs) - ] - ) - bounds = array( - [ - [self.constraint_lbs[cn], self.constraint_ubs[cn]] - for cn in self.constraints - ] - ) - A = concatenate((A, Av)) - bounds = concatenate((bounds, vbounds)) - else: - A = Av - bounds = vbounds - if A.shape[0] == 0: - A = None - else: - A = csc_matrix( - (A[:, 2], (A[:, 0].astype("int64"), A[:, 1].astype("int64"))), - shape=(nc + nv, nv), - ) - return P, q, A, bounds + return settings def solve(self): """Solve the OSQP problem.""" - settings = self.settings.copy() + settings = self.map_settings() P, q, A, bounds = self.build() solver = osqp.OSQP() if P is None: # see https://github.com/cvxgrp/cvxpy/issues/898 settings.update({"adaptive_rho": 0, "rho": 1.0, "alpha": 1.0}) solver.setup(P=P, q=q, A=A, l=bounds[:, 0], u=bounds[:, 1], **settings) # noqa - if self.__solution is not None: + if self._solution is not None: if self.still_valid(A, bounds): - solver.warm_start(x=self.__solution["x"], y=self.__solution["y"]) - solver.update_settings(rho=self.__solution["rho"]) + solver.warm_start(x=self._solution["x"], y=self._solution["y"]) + solver.update_settings(rho=self._solution["rho"]) solution = solver.solve() nc = len(self.constraints) nv = len(self.variables) @@ -217,36 +92,26 @@ def solve(self): if nc > 0: self.cprimals = dict(zip(self.constraints, A.dot(solution.x)[0:nc])) self.duals = dict(zip(self.constraints, solution.y[0:nc])) - if not isnan(solution.info.obj_val): + if not np.isnan(solution.info.obj_val): self.obj_value = solution.info.obj_val * self.direction self.status = solution.info.status else: self.status = "primal infeasible" self.obj_value = None self.info = solution.info - self.__solution = { + self._solution = { "x": solution.x, "y": solution.y, "rho": solution.info.rho_estimate, } - def reset(self, everything=False): - """Reset the public solver solution.""" - self.info = None - self.primals = {} - self.cprimals = {} - self.duals = {} - self.vduals = {} - if everything: - self.__solution = None - def still_valid(self, A, bounds): """Check if previous solutions is still feasible.""" - if len(self.__solution["x"]) != len(self.variables) or len( - self.__solution["y"] + if len(self._solution["x"]) != len(self.variables) or len( + self._solution["y"] ) != len(self.constraints): return False - c = A.dot(self.__solution["x"]) + c = A.dot(self._solution["x"]) ea = self.settings["eps_abs"] er = self.settings["eps_rel"] valid = np.all( @@ -255,705 +120,27 @@ def still_valid(self, A, bounds): ) return valid - def clean(self): - """Remove unused variables and constraints.""" - self.reset() - self.variable_lbs = { - k: v for k, v in six.iteritems(self.variable_lbs) if k in self.variables - } - self.variable_ubs = { - k: v for k, v in six.iteritems(self.variable_ubs) if k in self.variables - } - self.constraint_coefs = { - k: v - for k, v in six.iteritems(self.constraint_coefs) - if k[0] in self.constraints and k[1] in self.variables - } - self.obj_linear_coefs = { - k: v for k, v in six.iteritems(self.obj_linear_coefs) if k in self.variables - } - self.obj_quadratic_coefs = { - k: v - for k, v in six.iteritems(self.obj_quadratic_coefs) - if k[0] in self.variables and k[1] in self.variables - } - self.constraint_lbs = { - k: v for k, v in six.iteritems(self.constraint_lbs) if k in self.constraints - } - self.constraint_ubs = { - k: v for k, v in six.iteritems(self.constraint_ubs) if k in self.constraints - } - - def rename_constraint(self, old, new): - """Rename a constraint.""" - self.reset() - self.constraints.remove(old) - self.constraints.add(new) - name_map = {k: k for k in self.constraints} - name_map[old] = new - self.constraint_coefs = { - (name_map[k[0]], k[1]): v for k, v in six.iteritems(self.constraint_coefs) - } - self.constraint_lbs[new] = self.constraint_lbs.pop(old) - self.constraint_ubs[new] = self.constraint_ubs.pop(old) - - def rename_variable(self, old, new): - """Rename a variable.""" - self.reset() - self.variables.remove(old) - self.variables.add(new) - name_map = {k: k for k in self.variables} - name_map[old] = new - self.constraint_coefs = { - (k[0], name_map[k[1]]): v for k, v in six.iteritems(self.constraint_coefs) - } - self.obj_quadratic_coefs = { - (name_map[k[0]], name_map[k[1]]): v - for k, v in six.iteritems(self.obj_quadratic_coefs) - } - self.obj_linear_coefs = { - name_map[k]: v for k, v in six.iteritems(self.obj_linear_coefs) - } - self.variable_lbs[new] = self.variable_lbs.pop(old) - self.variable_ubs[new] = self.variable_ubs.pop(old) - @six.add_metaclass(inheritdocstring) -class Variable(interface.Variable): - def __init__(self, name, *args, **kwargs): - super(Variable, self).__init__(name, **kwargs) - - @interface.Variable.type.setter - def type(self, value): - if self.problem is not None: - if value not in _TYPES: - raise ValueError( - "OSQP cannot handle variables of type '%s'. " % value - + "The following variable types are available: " - + ", ".join(_TYPES) - ) - super(Variable, Variable).type.fset(self, value) - - def _get_primal(self): - return self.problem.problem.primals.get(self.name, None) - - @property - def dual(self): - if self.problem is not None: - return self.problem.problem.vduals.get(self.name, None) - return None - - @interface.Variable.name.setter - def name(self, value): - old_name = getattr(self, "name", None) - super(Variable, Variable).name.fset(self, value) - if getattr(self, "problem", None) is not None: - if old_name != value: - self.problem.problem.rename_variable(old_name, value) +class Variable(mi.Variable): + pass @six.add_metaclass(inheritdocstring) -class Constraint(interface.Constraint): +class Constraint(mi.Constraint): _INDICATOR_CONSTRAINT_SUPPORT = False - def __init__(self, expression, sloppy=False, *args, **kwargs): - super(Constraint, self).__init__(expression, *args, sloppy=sloppy, **kwargs) - - def set_linear_coefficients(self, coefficients): - if self.problem is not None: - self.problem.update() - self.problem.problem.reset() - for var, coef in six.iteritems(coefficients): - self.problem.problem.constraint_coefs[(self.name, var.name)] = float( - coef - ) - else: - raise Exception( - "Can't change coefficients if constraint is not associated " - "with a model." - ) - - def get_linear_coefficients(self, variables): - if self.problem is not None: - self.problem.update() - coefs = { - v: self.problem.problem.constraint_coefs.get((self.name, v.name), 0.0) - for v in variables - } - return coefs - else: - raise Exception( - "Can't get coefficients from solver if constraint is not " "in a model" - ) - - def _get_expression(self): - if self.problem is not None: - variables = self.problem._variables - all_coefs = self.problem.problem.constraint_coefs - coefs = [(v, all_coefs.get((self.name, v.name), 0.0)) for v in variables] - expression = add([mul((symbolics.Real(co), v)) for (v, co) in coefs]) - self._expression = expression - return self._expression - - @property - def problem(self): - return self._problem - - @problem.setter - def problem(self, value): - if value is None: - # Update expression from solver instance one last time - self._get_expression() - self._problem = value - - @property - def primal(self): - if self.problem is None: - return None - return self.problem.problem.cprimals.get(self.name, None) - - @property - def dual(self): - if self.problem is None: - return None - d = self.problem.problem.duals.get(self.name, None) - if d is not None: - d = -d - return d - - @interface.Constraint.name.setter - def name(self, value): - old_name = self.name - super(Constraint, Constraint).name.fset(self, value) - if getattr(self, "problem", None) is not None: - if old_name != value: - self.problem.problem.rename_constraint(old_name, value) - - @interface.Constraint.lb.setter - def lb(self, value): - self._check_valid_lower_bound(value) - if getattr(self, "problem", None) is not None: - lb = -Infinity if value is None else float(value) - self.problem.problem.constraint_lbs[self.name] = lb - self._lb = value - - @interface.Constraint.ub.setter - def ub(self, value): - self._check_valid_upper_bound(value) - if getattr(self, "problem", None) is not None: - ub = Infinity if value is None else float(value) - self.problem.problem.constraint_ubs[self.name] = ub - self._ub = value - - def __iadd__(self, other): - # if self.problem is not None: - # self.problem._add_to_constraint(self.index, other) - if self.problem is not None: - problem_reference = self.problem - self.problem._remove_constraint(self) - super(Constraint, self).__iadd__(other) - problem_reference._add_constraint(self, sloppy=False) - else: - super(Constraint, self).__iadd__(other) - return self - @six.add_metaclass(inheritdocstring) -class Objective(interface.Objective): - def __init__(self, expression, sloppy=False, **kwargs): - super(Objective, self).__init__(expression, sloppy=sloppy, **kwargs) - self._expression_expired = False - if not (sloppy or self.is_Linear or self.is_Quadratic): - raise ValueError("OSQP only supports linear and quadratic objectives.") - - @property - def value(self): - if getattr(self, "problem", None) is None: - return None - if self.problem.problem.obj_value is None: - return None - return self.problem.problem.obj_value + self.problem._objective_offset - - @interface.Objective.direction.setter - def direction(self, value): - super(Objective, Objective).direction.__set__(self, value) - if value == "min": - self.problem.problem.direction = 1 - else: - self.problem.problem.direction = -1 - - def _get_expression(self): - if ( - self.problem is not None - and self._expression_expired - and len(self.problem._variables) > 0 - ): - model = self.problem - vars = model._variables - expression = add( - [ - coef * vars[vn] - for vn, coef in six.iteritems(model.problem.obj_linear_coefs) - ] - ) - q_ex = add( - [ - coef * vars[vn[0]] * vars[vn[1]] - for vn, coef in six.iteritems(model.problem.obj_quadratic_coefs) - ] - ) - expression += q_ex - self._expression = expression - self._expression_expired = False - return self._expression - - def set_linear_coefficients(self, coefficients): - if self.problem is not None: - self.problem.update() - for v, coef in six.iteritems(coefficients): - self.problem.problem.obj_linear_coefs[v.name] = float(coef) - self._expression_expired = True - else: - raise Exception( - "Can't change coefficients if objective is not associated " - "with a model." - ) - - def get_linear_coefficients(self, variables): - if self.problem is not None: - self.problem.update() - coefs = { - v: self.problem.problem.obj_linear_coefs.get(v.name, 0.0) - for v in variables - } - return coefs - else: - raise Exception( - "Can't get coefficients from solver if objective " "is not in a model" - ) +class Objective(mi.Objective): + pass @six.add_metaclass(inheritdocstring) -class Configuration(interface.MathematicalProgrammingConfiguration): - def __init__( - self, - lp_method="primal", - presolve=False, - verbosity=0, - timeout=None, - qp_method="primal", - linear_solver="qdldl", - *args, - **kwargs - ): - super(Configuration, self).__init__(*args, **kwargs) - self.lp_method = lp_method - self.presolve = presolve - self.verbosity = verbosity - self.timeout = timeout - self.qp_method = qp_method - self.linear_solver = linear_solver - if "tolerances" in kwargs: - for key, val in six.iteritems(kwargs["tolerances"]): - if key in self._tolerance_functions(): - setattr(self.tolerances, key, val) - - @property - def lp_method(self): - """The algorithm used to solve LP problems.""" - return "primal" - - @lp_method.setter - def lp_method(self, lp_method): - if lp_method not in _LP_METHODS: - raise ValueError( - "LP Method %s is not valid (choose one of: %s)" - % (lp_method, ", ".join(_LP_METHODS)) - ) - - @property - def linear_solver(self): - return self._linear_solver - - @linear_solver.setter - def linear_solver(self, solver): - if solver not in ("qdldl", "mkl pardiso"): - raise ValueError( - "%s is not valid (choose either `qdldl` or `mkl pardiso`)" % solver - ) - if getattr(self, "problem", None) is not None: - self.problem.problem.settings["linsys_solver"] = solver - self._linear_solver = solver - - def _set_presolve(self, value): - if getattr(self, "problem", None) is not None: - if value is True: - self.problem.problem.settings["scaling"] = 10 - elif value is False or value == "auto": - self.problem.problem.settings["scaling"] = 0 - else: - raise ValueError( - str(value) + " is not a valid presolve parameter. Must be True, " - "False or 'auto'." - ) - - @property - def presolve(self): - return self.problem.problem.settings["scaling"] > 0 - - @presolve.setter - def presolve(self, value): - self._set_presolve(value) - self._presolve = value - - @property - def verbosity(self): - return self._verbosity - - @verbosity.setter - def verbosity(self, value): - if getattr(self, "problem", None) is not None: - self.problem.problem.settings["verbose"] = int(value > 0) - self._verbosity = value - - @property - def timeout(self): - return self.problem.problem.settings["time_limit"] - - @timeout.setter - def timeout(self, value): - if getattr(self, "problem", None) is not None: - if value is None: - self.problem.problem.settings["time_limit"] = 0 - else: - self.problem.problem.settings["time_limit"] = value - - def __getstate__(self): - return { - "presolve": self.presolve, - "timeout": self.timeout, - "verbosity": self.verbosity, - "linear_solver": self.linear_solver, - "tolerances": { - "feasibility": self.tolerances.feasibility, - "optimality": self.tolerances.optimality, - "integrality": self.tolerances.integrality, - }, - } - - def __setstate__(self, state): - for key, val in six.iteritems(state): - if key != "tolerances": - setattr(self, key, val) - for key, val in six.iteritems(state["tolerances"]): - if key in self._tolerance_functions(): - setattr(self.tolerances, key, val) - - @property - def qp_method(self): - """Change the algorithm used to optimize QP problems.""" - return "primal" - - @qp_method.setter - def qp_method(self, value): - if value not in _QP_METHODS: - raise ValueError( - "%s is not a valid qp_method. Choose between %s" - % (value, str(_QP_METHODS)) - ) - self._qp_method = value - - def _get_feasibility(self): - return self.problem.problem.settings["eps_prim_inf"] - - def _set_feasibility(self, value): - self.problem.problem.settings["eps_prim_inf"] = value - self.problem.problem.settings["eps_dual_inf"] = value - - def _get_integrality(self): - return 1e-6 - - def _set_integrality(self, value): - pass - - def _get_optimality(self): - return self.problem.problem.settings["eps_abs"] - - def _set_optimality(self, value): - self.problem.problem.settings["eps_abs"] = value - self.problem.problem.settings["eps_rel"] = value - - def _tolerance_functions(self): - return { - "feasibility": (self._get_feasibility, self._set_feasibility), - "optimality": (self._get_optimality, self._set_optimality), - "integrality": (self._get_integrality, self._set_integrality), - } +class Configuration(mi.Configuration): + pass @six.add_metaclass(inheritdocstring) -class Model(interface.Model): - def _initialize_problem(self): - self.problem = OSQPProblem() - - def _initialize_model_from_problem(self, problem, vc_mapping=None, offset=0): - if not isinstance(problem, OSQPProblem): - raise TypeError("Provided problem is not a valid OSQP model.") - self.problem = problem - for name in self.problem.variables: - var = Variable( - name, - lb=self.problem.variable_lbs[name], - ub=self.problem.variable_ubs[name], - problem=self, - ) - super(Model, self)._add_variables([var]) - - for name in self.problem.constraints: - # Since constraint expressions are lazily retrieved from the - # solver they don't have to be built here - lhs = symbolics.Integer(0) - constr = Constraint( - lhs, - lb=self.problem.constraint_lbs[name], - ub=self.problem.constraint_ubs[name], - name=name, - problem=self, - ) - - super(Model, self)._add_constraints([constr], sloppy=True) - - if vc_mapping is None: - for constr in self.constraints: - name = constr.name - for variable in constr.variables: - try: - self._variables_to_constraints_mapping[variable.name].add(name) - except KeyError: - self._variables_to_constraints_mapping[variable.name] = set( - [name] - ) - else: - self._variables_to_constraints_mapping = vc_mapping - - linear_expression = add( - [ - coef * self._variables[vn] - for vn, coef in six.iteritems(self.problem.obj_linear_coefs) - ] - ) - quadratic_expression = add( - [ - coef * self._variables[vn[0]] * self._variables[vn[1]] - for vn, coef in six.iteritems(self.problem.obj_quadratic_coefs) - ] - ) - - self._objective_offset = offset - self._objective = Objective( - linear_expression + quadratic_expression + offset, - problem=self, - direction={-1: "max", 1: "min"}[self.problem.direction], - name="osqp_objective", - ) - - @property - def objective(self): - return self._objective - - @objective.setter - def objective(self, value): - value.problem = None - if self._objective is not None: # Reset previous objective - self.problem.obj_linear_coefs = dict() - self.problem.obj_quadratic_coefs = dict() - super(Model, self.__class__).objective.fset(self, value) - self.update() - expression = self._objective._expression - ( - offset, - linear_coefficients, - quadratic_coeffients, - ) = parse_optimization_expression(value, quadratic=True, expression=expression) - self._objective_offset = offset - if linear_coefficients: - self.problem.obj_linear_coefs = { - v.name: float(c) for v, c in six.iteritems(linear_coefficients) - } - - for key, coef in six.iteritems(quadratic_coeffients): - if len(key) == 1: - var = six.next(iter(key)) - self.problem.obj_quadratic_coefs[(var.name, var.name)] = float(coef) - else: - var1, var2 = key - self.problem.obj_quadratic_coefs[(var1.name, var2.name)] = 0.5 * float( - coef - ) - self.problem.obj_quadratic_coefs[(var2.name, var1.name)] = 0.5 * float( - coef - ) - - self._set_objective_direction(value.direction) - value.problem = self - - def _set_objective_direction(self, direction): - self.problem.direction = {"min": 1, "max": -1}[direction] - - def _get_primal_values(self): - if len(self.problem.primals) == 0: - raise SolverError("The problem has not been solved yet!") - primal_values = [self.problem.primals[v.name] for v in self._variables] - return primal_values - - def _get_reduced_costs(self): - if len(self.problem.duals) == 0: - raise SolverError("The problem has not been solved yet!") - reduced_costs = [self.problem.vduals[v.name] for v in self._variables] - return reduced_costs - - def _get_constraint_values(self): - if len(self.problem.primals) == 0: - raise SolverError("The problem has not been solved yet!") - constraint_primals = [self.problem.cprimals[c.name] for c in self._constraints] - return constraint_primals - - def _get_shadow_prices(self): - if len(self.problem.primals) == 0: - raise SolverError("The problem has not been solved yet!") - dual_values = [-self.problem.duals[c.name] for c in self._constraints] - return dual_values - - @property - def is_integer(self): - return False - - def _optimize(self): - self.update() - self.problem.solve() - osqp_status = self.problem.status - self._original_status = osqp_status - status = _STATUS_MAP[osqp_status] - return status - - def _set_variable_bounds_on_problem(self, var_lb, var_ub): - self.problem.reset() - for var, val in var_lb: - lb = -Infinity if val is None else float(val) - self.problem.variable_lbs[var.name] = lb - for var, val in var_ub: - ub = Infinity if val is None else val - self.problem.variable_ubs[var.name] = float(ub) - - def _add_variables(self, variables): - super(Model, self)._add_variables(variables) - self.problem.reset() - for variable in variables: - lb = -Infinity if variable.lb is None else float(variable.lb) - ub = Infinity if variable.ub is None else float(variable.ub) - self.problem.variables.add(variable.name) - self.problem.variable_lbs[variable.name] = lb - self.problem.variable_ubs[variable.name] = ub - variable.problem = self - - def _remove_variables(self, variables): - # Not calling parent method to avoid expensive variable removal from sympy expressions - if self.objective is not None: - self.objective._expression = self.objective.expression.xreplace( - {var: 0 for var in variables} - ) - for variable in variables: - del self._variables_to_constraints_mapping[variable.name] - variable.problem = None - del self._variables[variable.name] - self.problem.variables.remove(variable.name) - self.problem.clean() - - def _add_constraints(self, constraints, sloppy=False): - super(Model, self)._add_constraints(constraints, sloppy=sloppy) - self.problem.reset() - for constraint in constraints: - constraint._problem = None # This needs to be done in order to not trigger constraint._get_expression() - if constraint.is_Linear: - _, coeff_dict, _ = parse_optimization_expression(constraint) - lb = -Infinity if constraint.lb is None else float(constraint.lb) - ub = Infinity if constraint.ub is None else float(constraint.ub) - self.problem.constraints.add(constraint.name) - self.problem.constraint_coefs.update( - { - (constraint.name, v.name): float(co) - for v, co in six.iteritems(coeff_dict) - } - ) - self.problem.constraint_lbs[constraint.name] = lb - self.problem.constraint_ubs[constraint.name] = ub - constraint.problem = self - elif constraint.is_Quadratic: - raise NotImplementedError( - "Quadratic constraints (like %s) are not supported " - "in OSQP yet." % constraint - ) - else: - raise ValueError( - "OSQP only supports linear or quadratic constraints. " - "%s is neither linear nor quadratic." % constraint - ) - - def _remove_constraints(self, constraints): - super(Model, self)._remove_constraints(constraints) - for constraint in constraints: - self.problem.constraints.remove(constraint.name) - self.problem.clean() - - def _get_variable_indices(self, names): - vmap = dict(zip(self.variables, range(len(self.variables)))) - return [vmap[n] for n in names] - - def __setstate__(self, d): - self.__init__() - osqp = pickle.loads(d["osqp"]) - # workaround to conserve the order - osqp.variables = d["variables"] - osqp.constraints = d["constraints"] - self._initialize_model_from_problem( - osqp, vc_mapping=d["v_to_c"], offset=d["offset"] - ) - osqp.variables = set(osqp.variables) - osqp.constraints = set(osqp.constraints) - self.configuration = Configuration() - self.configuration.problem = self - self.configuration.__setstate__(d["config"]) - - def __getstate__(self): - self.problem.reset() - self.update() - return { - "osqp": pickle.dumps(self.problem), - "variables": [v.name for v in self._variables], - "constraints": [c.name for c in self._constraints], - "v_to_c": self._variables_to_constraints_mapping, - "config": self.configuration.__getstate__(), - "offset": getattr(self, "_objective_offset", 0.0), - } - - @classmethod - def from_lp(self, lp_problem_str): - """Read a model from an LP file. - - OSQP does not have an integrated LP reader so it will either use - cplex or glpk to read the model. This means that QP problems will - currently require cplex to be read :( - """ - if available_solvers["CPLEX"]: - from optlang import cplex_interface - - mod = cplex_interface.Model.from_lp(lp_problem_str) - mod.configuration.lp_method = "auto" - mod.configuration.qp_method = "auto" - return super(Model, self).clone(mod) - else: - from optlang import glpk_interface - - mod = glpk_interface.Model.from_lp(lp_problem_str) - mod.configuration.lp_method = "auto" - return super(Model, self).clone(mod) +class Model(mi.Model): + ProblemClass = OSQPProblem diff --git a/src/optlang/tests/test_osqp_interface.py b/src/optlang/tests/test_osqp_interface.py index 6cb5e43..7f7aa5e 100644 --- a/src/optlang/tests/test_osqp_interface.py +++ b/src/optlang/tests/test_osqp_interface.py @@ -475,25 +475,17 @@ def test_tolerance_parameters(self): getattr(model.configuration.tolerances, param), 2 * val ) - def test_solver(self): - for option in ("qdldl", "mkl pardiso"): - self.configuration.linear_solver = option - self.assertEqual(self.configuration.linear_solver, option) - self.assertEqual(self.model.problem.settings["linsys_solver"], option) - - self.assertRaises(ValueError, setattr, self.configuration, "lp_method", "weird_stuff") - def test_lp_method(self): - for option in ("auto", "primal"): + for option in ("auto", ): self.configuration.lp_method = option - self.assertEqual(self.configuration.lp_method, "primal") + self.assertEqual(self.configuration.lp_method, "auto") self.assertRaises(ValueError, setattr, self.configuration, "lp_method", "weird_stuff") def test_qp_method(self): - for option in osqp_interface._QP_METHODS: + for option in ("auto", ): self.configuration.qp_method = option - self.assertEqual(self.configuration.qp_method, "primal") + self.assertEqual(self.configuration.qp_method, "auto") self.assertRaises(ValueError, setattr, self.configuration, "qp_method", "weird_stuff") From 56ae461b0e615e9fb180cbbb096bb6aeb244fb36 Mon Sep 17 00:00:00 2001 From: Christian Diener Date: Fri, 1 Sep 2023 11:56:24 -0700 Subject: [PATCH 04/10] some cleanup --- src/optlang/matrix_interface.py | 177 ++++++++++++++++++++++++-------- src/optlang/osqp_interface.py | 28 ++++- 2 files changed, 159 insertions(+), 46 deletions(-) diff --git a/src/optlang/matrix_interface.py b/src/optlang/matrix_interface.py index 1ab2449..e5cb4bc 100644 --- a/src/optlang/matrix_interface.py +++ b/src/optlang/matrix_interface.py @@ -15,12 +15,26 @@ This interface is only supposed to be a basis to add matrix solvers to optlang. It creates a thin layer between all optlang operations and the standard formulation of LP, QP, and MILP problems. + +To create a solver from this, you need to do the following: + +1. Create a class that inherits from MatrixProblem and implement the `solve` method + and optionally overwrite the `reset` and `still_valid` methods. + +2. Create derived classes for Variable, Constraint, Objective, and Configuration. + You don't need to overwrite anything. Those should work out of the box. + +3. Create a derived class for Model and set the `ProblemClass` attribute to the class + from (1) and the `status_map` to map the solver status to an Optlang status from + `optlang.interface`. + +For an example take a look at `optlang.osqp_interface`. """ import logging import six import pickle -import numpy as np +from collections import defaultdict from numpy import array, concatenate, Infinity, zeros from optlang import interface, symbolics, available_solvers @@ -35,21 +49,7 @@ log = logging.getLogger(__name__) -_STATUS_MAP = { - "interrupted by user": interface.ABORTED, - "run time limit reached": interface.TIME_LIMIT, - "feasible": interface.FEASIBLE, - "primal infeasible": interface.INFEASIBLE, - "dual infeasible": interface.INFEASIBLE, - "primal infeasible inaccurate": interface.INFEASIBLE, - "dual infeasible inaccurate": interface.INFEASIBLE, - "solved inaccurate": interface.NUMERIC, - "solved": interface.OPTIMAL, - "maximum iterations reached": interface.ITERATION_LIMIT, - "unsolved": interface.SPECIAL, - "problem non convex": interface.SPECIAL, - "non-existing-status": "Here for testing that missing statuses are handled.", -} +_STATUS_MAP = defaultdict(lambda: interface.UNDEFINED) _LP_METHODS = ("auto", ) @@ -63,8 +63,7 @@ class MatrixProblem(object): This class assumes that the problem will be pretty much immutable. This is a small intermediate layer based on dictionaries that is fast to modify - but can also be converted to a matrix formulation quickly. This is the only - thing you will need to modify. + but can also be converted to a matrix formulation quickly. """ def __init__(self): @@ -90,7 +89,13 @@ def __init__(self): self.default_settings() def default_settings(self): - """Set the default settings of the solver.""" + """Set the default settings of the solver. + + Returns + ------- + dict + The default settings for the solver. + """ self.settings = { "presolve": True, "lp_method": "auto", @@ -99,11 +104,33 @@ def default_settings(self): "primal_inf_tolerance": 1e-6, "dual_inf_tolerance": 1e-6, "optimality_tolerance": 1e-6, - "verbose": 0 + "mip_tolerance": 1e-6, + "verbose": 1, + "threads": 1, + "iteration_limit": 100000 } - def build(self): - """Build the problem instance.""" + def build(self, add_variable_constraints=False): + """Build the problem instance. + + Parameters + ---------- + add_variable_constraints : bool + Whether to add constraints for each individual variables as a (sparse) + block diagonal to the constraint matrix A. + + Returns + ------- + tuple of [P, q, A, bounds, vbounds, integer] with the following components: + - P : Quadratic objective coefficients as an upper triangular CSC matrix + - q : Linear objective coefficients as an array. + - A : Constraint matrix as a CSC matrix. + - bounds : Constraint lower and upper bounds as a 2xNc matrix. + - vbounds : Variable bounds as a 2xNv matrix. + - integer : Binary array indicating whether variable i is integer + (for MIPs). + + """ d = self.direction vmap = dict(zip(self.variables, range(len(self.variables)))) nv = len(self.variables) @@ -127,10 +154,10 @@ def build(self): self.obj_linear_coefs.values() ) q = q * d - Av = array([[vmap[k] + nc, vmap[k], 1.0] for k in self.variables]) vbounds = array( [[self.variable_lbs[vn], self.variable_ubs[vn]] for vn in self.variables] ) + if len(self.constraint_coefs) > 0: A = array( [ @@ -144,26 +171,51 @@ def build(self): for cn in self.constraints ] ) - A = concatenate((A, Av)) - bounds = concatenate((bounds, vbounds)) - else: - A = Av - bounds = vbounds - if A.shape[0] == 0: - A = None - else: + if add_variable_constraints: + Av = array([[vmap[k] + nc, vmap[k], 1.0] for k in self.variables]) + A = concatenate((A, Av)) + bounds = concatenate((bounds, vbounds)) A = csc_matrix( (A[:, 2], (A[:, 0].astype("int64"), A[:, 1].astype("int64"))), shape=(nc + nv, nv), ) - return P, q, A, bounds + elif add_variable_constraints: + Av = array([[vmap[k], vmap[k], 1.0] for k in self.variables]) + A = csc_matrix( + (Av[:, 2], (Av[:, 0].astype("int64"), Av[:, 1].astype("int64"))), + shape=(nv, nv), + ) + bounds = vbounds + else: + A = None + bounds = None + + integer = array([int(vn in self.integer_vars) for vn in self.variables]) + return P, q, A, bounds, vbounds, integer def solve(self): - """Solve the problem.""" + """Solve the problem. + + This will return nothing but will fill the solution attributes. + + Returns + ------- + Nothing. + """ raise NotImplementedError("This needs to be overwritten by the child class.") def reset(self, everything=False): - """Reset the public solver solution.""" + """Reset the public solver solution. + + Parameters + ---------- + everything : bool + If true will also clear cached solution bases or warm start solutions. + + Returns + ------- + Nothing. + """ self.info = None self.primals = {} self.cprimals = {} @@ -173,7 +225,21 @@ def reset(self, everything=False): self._solution = None def still_valid(self, A, bounds): - """Check if previous solutions is still feasible.""" + """Check if previous solutions is still feasible. + + Parameters + ---------- + A : numpy.csc_matrix + The constraint matrix A. + bounds : numpy.array + The constraint bounds. + + Returns + ------- + bool + Whether the current solution still fulfills the bounds. + + """ raise NotImplementedError("This needs to be overwritten by the child class.") def clean(self): @@ -207,7 +273,20 @@ def clean(self): self.integer_vars = {v for v in self.integer_vars if v in self.variables} def rename_constraint(self, old, new): - """Rename a constraint.""" + """Rename a constraint. + + Parameters + ---------- + old : str + The old name of the constraint. + new : str + The new name of the constraint. + + Returns + ------- + Nothing. + + """ self.reset() self.constraints.remove(old) self.constraints.add(new) @@ -220,7 +299,20 @@ def rename_constraint(self, old, new): self.constraint_ubs[new] = self.constraint_ubs.pop(old) def rename_variable(self, old, new): - """Rename a variable.""" + """Rename a variable. + + Parameters + ---------- + old : str + The old name of the variable. + new : str + The new name of the variable. + + Returns + ------- + Nothing. + + """ self.reset() self.variables.remove(old) self.variables.add(new) @@ -494,10 +586,8 @@ def lp_method(self, lp_method): def _set_presolve(self, value): if getattr(self, "problem", None) is not None: - if isinstance(value, bool): + if value in [True, False, "auto"]: self.problem.problem.settings["presolve"] = value - elif value == "auto": - self.problem.problem.settings["presove"] = False else: raise ValueError( str(value) + " is not a valid presolve parameter. Must be True, " @@ -599,6 +689,7 @@ def _tolerance_functions(self): @six.add_metaclass(inheritdocstring) class Model(interface.Model): ProblemClass = MatrixProblem + status_map = _STATUS_MAP def _initialize_problem(self): self.problem = self.ProblemClass() @@ -733,14 +824,14 @@ def _get_shadow_prices(self): @property def is_integer(self): - return False + return len(self.problem.integer_vars) > 0 def _optimize(self): self.update() self.problem.solve() prior_status = self.problem.status self._original_status = prior_status - status = _STATUS_MAP[prior_status] + status = self.status_map[prior_status] return status def _set_variable_bounds_on_problem(self, var_lb, var_ub): @@ -761,6 +852,8 @@ def _add_variables(self, variables): self.problem.variables.add(variable.name) self.problem.variable_lbs[variable.name] = lb self.problem.variable_ubs[variable.name] = ub + if variable.type in ("binary", "integer"): + self.problem.integer_vars.add(variable.name) variable.problem = self def _remove_variables(self, variables): diff --git a/src/optlang/osqp_interface.py b/src/optlang/osqp_interface.py index 4e68033..3bfa5d0 100644 --- a/src/optlang/osqp_interface.py +++ b/src/optlang/osqp_interface.py @@ -24,6 +24,7 @@ import numpy as np import optlang.matrix_interface as mi +from optlang import interface from optlang.util import inheritdocstring log = logging.getLogger(__name__) @@ -43,6 +44,23 @@ raise ImportError("The osqp_interface requires osqp or cuosqp!") +_STATUS_MAP = { + "interrupted by user": interface.ABORTED, + "run time limit reached": interface.TIME_LIMIT, + "feasible": interface.FEASIBLE, + "primal infeasible": interface.INFEASIBLE, + "dual infeasible": interface.INFEASIBLE, + "primal infeasible inaccurate": interface.INFEASIBLE, + "dual infeasible inaccurate": interface.INFEASIBLE, + "solved inaccurate": interface.NUMERIC, + "solved": interface.OPTIMAL, + "maximum iterations reached": interface.ITERATION_LIMIT, + "unsolved": interface.SPECIAL, + "problem non convex": interface.SPECIAL, + "non-existing-status": "Here for testing that missing statuses are handled." +} + + class OSQPProblem(mi.MatrixProblem): """A concise representation of an OSQP problem. @@ -51,7 +69,7 @@ class OSQPProblem(mi.MatrixProblem): but can also be converted to an OSQP problem without too much hassle. """ - def map_settings(self): + def osqp_settings(self): """Map internal settings to OSQP settings.""" settings = { "linsys_solver": "qdldl", @@ -72,8 +90,8 @@ def map_settings(self): def solve(self): """Solve the OSQP problem.""" - settings = self.map_settings() - P, q, A, bounds = self.build() + settings = self.osqp_settings() + P, q, A, bounds, _, _ = self.build(add_variable_constraints=True) solver = osqp.OSQP() if P is None: # see https://github.com/cvxgrp/cvxpy/issues/898 @@ -82,7 +100,8 @@ def solve(self): if self._solution is not None: if self.still_valid(A, bounds): solver.warm_start(x=self._solution["x"], y=self._solution["y"]) - solver.update_settings(rho=self._solution["rho"]) + if "rho" in self._solution: + solver.update_settings(rho=self._solution["rho"]) solution = solver.solve() nc = len(self.constraints) nv = len(self.variables) @@ -144,3 +163,4 @@ class Configuration(mi.Configuration): @six.add_metaclass(inheritdocstring) class Model(mi.Model): ProblemClass = OSQPProblem + status_map = _STATUS_MAP From 53d77e2c493ce1cc32977998a40bfb60683b389f Mon Sep 17 00:00:00 2001 From: Christian Diener Date: Mon, 11 Sep 2023 13:46:37 -0700 Subject: [PATCH 05/10] review comments --- src/optlang/matrix_interface.py | 85 +++++++++++++++++---------------- src/optlang/osqp_interface.py | 22 +++++---- 2 files changed, 56 insertions(+), 51 deletions(-) diff --git a/src/optlang/matrix_interface.py b/src/optlang/matrix_interface.py index e5cb4bc..e0bd6d7 100644 --- a/src/optlang/matrix_interface.py +++ b/src/optlang/matrix_interface.py @@ -31,14 +31,13 @@ For an example take a look at `optlang.osqp_interface`. """ +import abc import logging -import six import pickle -from collections import defaultdict +from collections import defaultdict, namedtuple from numpy import array, concatenate, Infinity, zeros from optlang import interface, symbolics, available_solvers -from optlang.util import inheritdocstring from optlang.expression_parsing import parse_optimization_expression from optlang.exceptions import SolverError @@ -57,8 +56,14 @@ _TYPES = ("continuous", "binary", "integer") +SparseProblem = namedtuple( + "SparseProblem", + ["P", "q", "A", "bounds", "vbounds", "integer"] +) +"""A representation of the Problem in standard form.""" -class MatrixProblem(object): + +class MatrixProblem(abc.ABC): """A concise representation of an LP/QP/MILP problem. This class assumes that the problem will be pretty much immutable. This is @@ -66,7 +71,8 @@ class MatrixProblem(object): but can also be converted to a matrix formulation quickly. """ - def __init__(self): + def __init__(self, **kwargs): + super().__init__(**kwargs) self.variables = set() self.integer_vars = set() self.constraints = set() @@ -140,7 +146,7 @@ def build(self, add_variable_constraints=False): P = array( [ [vmap[vn[0]], vmap[vn[1]], coef * d * 2.0] - for vn, coef in six.iteritems(self.obj_quadratic_coefs) + for vn, coef in self.obj_quadratic_coefs.items() ] ) P = csc_matrix( @@ -162,7 +168,7 @@ def build(self, add_variable_constraints=False): A = array( [ [cmap[vn[0]], vmap[vn[1]], coef] - for vn, coef in six.iteritems(self.constraint_coefs) + for vn, coef in self.constraint_coefs.items() ] ) bounds = array( @@ -191,8 +197,9 @@ def build(self, add_variable_constraints=False): bounds = None integer = array([int(vn in self.integer_vars) for vn in self.variables]) - return P, q, A, bounds, vbounds, integer + return SparseProblem(P, q, A, bounds, vbounds, integer) + @abc.abstractmethod def solve(self): """Solve the problem. @@ -224,6 +231,7 @@ def reset(self, everything=False): if everything: self._solution = None + @abc.abstractmethod def still_valid(self, A, bounds): """Check if previous solutions is still feasible. @@ -242,33 +250,33 @@ def still_valid(self, A, bounds): """ raise NotImplementedError("This needs to be overwritten by the child class.") - def clean(self): + def prune(self): """Remove unused variables and constraints.""" self.reset() self.variable_lbs = { - k: v for k, v in six.iteritems(self.variable_lbs) if k in self.variables + k: v for k, v in self.variable_lbs.items() if k in self.variables } self.variable_ubs = { - k: v for k, v in six.iteritems(self.variable_ubs) if k in self.variables + k: v for k, v in self.variable_ubs.items() if k in self.variables } self.constraint_coefs = { k: v - for k, v in six.iteritems(self.constraint_coefs) + for k, v in self.constraint_coefs.items() if k[0] in self.constraints and k[1] in self.variables } self.obj_linear_coefs = { - k: v for k, v in six.iteritems(self.obj_linear_coefs) if k in self.variables + k: v for k, v in self.obj_linear_coefs.items() if k in self.variables } self.obj_quadratic_coefs = { k: v - for k, v in six.iteritems(self.obj_quadratic_coefs) + for k, v in self.obj_quadratic_coefs.items() if k[0] in self.variables and k[1] in self.variables } self.constraint_lbs = { - k: v for k, v in six.iteritems(self.constraint_lbs) if k in self.constraints + k: v for k, v in self.constraint_lbs.items() if k in self.constraints } self.constraint_ubs = { - k: v for k, v in six.iteritems(self.constraint_ubs) if k in self.constraints + k: v for k, v in self.constraint_ubs.items() if k in self.constraints } self.integer_vars = {v for v in self.integer_vars if v in self.variables} @@ -293,7 +301,7 @@ def rename_constraint(self, old, new): name_map = {k: k for k in self.constraints} name_map[old] = new self.constraint_coefs = { - (name_map[k[0]], k[1]): v for k, v in six.iteritems(self.constraint_coefs) + (name_map[k[0]], k[1]): v for k, v in self.constraint_coefs.items() } self.constraint_lbs[new] = self.constraint_lbs.pop(old) self.constraint_ubs[new] = self.constraint_ubs.pop(old) @@ -319,23 +327,22 @@ def rename_variable(self, old, new): name_map = {k: k for k in self.variables} name_map[old] = new self.constraint_coefs = { - (k[0], name_map[k[1]]): v for k, v in six.iteritems(self.constraint_coefs) + (k[0], name_map[k[1]]): v for k, v in self.constraint_coefs.items() } self.obj_quadratic_coefs = { (name_map[k[0]], name_map[k[1]]): v - for k, v in six.iteritems(self.obj_quadratic_coefs) + for k, v in self.obj_quadratic_coefs.items() } self.obj_linear_coefs = { - name_map[k]: v for k, v in six.iteritems(self.obj_linear_coefs) + name_map[k]: v for k, v in self.obj_linear_coefs.items() } self.variable_lbs[new] = self.variable_lbs.pop(old) self.variable_ubs[new] = self.variable_ubs.pop(old) -@six.add_metaclass(inheritdocstring) class Variable(interface.Variable): def __init__(self, name, *args, **kwargs): - super(Variable, self).__init__(name, **kwargs) + super(Variable, self).__init__(name=name, **kwargs) @interface.Variable.type.setter def type(self, value): @@ -366,7 +373,6 @@ def name(self, value): self.problem.problem.rename_variable(old_name, value) -@six.add_metaclass(inheritdocstring) class Constraint(interface.Constraint): _INDICATOR_CONSTRAINT_SUPPORT = False @@ -377,7 +383,7 @@ def set_linear_coefficients(self, coefficients): if self.problem is not None: self.problem.update() self.problem.problem.reset() - for var, coef in six.iteritems(coefficients): + for var, coef in coefficients.items(): self.problem.problem.constraint_coefs[(self.name, var.name)] = float( coef ) @@ -472,7 +478,6 @@ def __iadd__(self, other): return self -@six.add_metaclass(inheritdocstring) class Objective(interface.Objective): def __init__(self, expression, sloppy=False, **kwargs): super(Objective, self).__init__(expression, sloppy=sloppy, **kwargs) @@ -508,13 +513,13 @@ def _get_expression(self): expression = add( [ coef * vars[vn] - for vn, coef in six.iteritems(model.problem.obj_linear_coefs) + for vn, coef in model.problem.obj_linear_coefs.items() ] ) q_ex = add( [ coef * vars[vn[0]] * vars[vn[1]] - for vn, coef in six.iteritems(model.problem.obj_quadratic_coefs) + for vn, coef in model.problem.obj_quadratic_coefs.items() ] ) expression += q_ex @@ -525,7 +530,7 @@ def _get_expression(self): def set_linear_coefficients(self, coefficients): if self.problem is not None: self.problem.update() - for v, coef in six.iteritems(coefficients): + for v, coef in coefficients.items(): self.problem.problem.obj_linear_coefs[v.name] = float(coef) self._expression_expired = True else: @@ -548,7 +553,6 @@ def get_linear_coefficients(self, variables): ) -@six.add_metaclass(inheritdocstring) class Configuration(interface.MathematicalProgrammingConfiguration): def __init__( self, @@ -567,7 +571,7 @@ def __init__( self.timeout = timeout self.qp_method = qp_method if "tolerances" in kwargs: - for key, val in six.iteritems(kwargs["tolerances"]): + for key, val in kwargs["tolerances"].items(): if key in self._tolerance_functions(): setattr(self.tolerances, key, val) @@ -638,10 +642,10 @@ def __getstate__(self): } def __setstate__(self, state): - for key, val in six.iteritems(state): + for key, val in state.items(): if key != "tolerances": setattr(self, key, val) - for key, val in six.iteritems(state["tolerances"]): + for key, val in state["tolerances"].items(): if key in self._tolerance_functions(): setattr(self.tolerances, key, val) @@ -686,7 +690,6 @@ def _tolerance_functions(self): } -@six.add_metaclass(inheritdocstring) class Model(interface.Model): ProblemClass = MatrixProblem status_map = _STATUS_MAP @@ -737,13 +740,13 @@ def _initialize_model_from_problem(self, problem, vc_mapping=None, offset=0): linear_expression = add( [ coef * self._variables[vn] - for vn, coef in six.iteritems(self.problem.obj_linear_coefs) + for vn, coef in self.problem.obj_linear_coefs.items() ] ) quadratic_expression = add( [ coef * self._variables[vn[0]] * self._variables[vn[1]] - for vn, coef in six.iteritems(self.problem.obj_quadratic_coefs) + for vn, coef in self.problem.obj_quadratic_coefs.items() ] ) @@ -776,12 +779,12 @@ def objective(self, value): self._objective_offset = offset if linear_coefficients: self.problem.obj_linear_coefs = { - v.name: float(c) for v, c in six.iteritems(linear_coefficients) + v.name: float(c) for v, c in linear_coefficients.items() } - for key, coef in six.iteritems(quadratic_coeffients): + for key, coef in quadratic_coeffients.items(): if len(key) == 1: - var = six.next(iter(key)) + var = next(iter(key)) self.problem.obj_quadratic_coefs[(var.name, var.name)] = float(coef) else: var1, var2 = key @@ -868,7 +871,7 @@ def _remove_variables(self, variables): variable.problem = None del self._variables[variable.name] self.problem.variables.remove(variable.name) - self.problem.clean() + self.problem.prune() def _add_constraints(self, constraints, sloppy=False): super(Model, self)._add_constraints(constraints, sloppy=sloppy) @@ -884,7 +887,7 @@ def _add_constraints(self, constraints, sloppy=False): self.problem.constraint_coefs.update( { (constraint.name, v.name): float(co) - for v, co in six.iteritems(coeff_dict) + for v, co in coeff_dict.items() } ) self.problem.constraint_lbs[constraint.name] = lb @@ -905,7 +908,7 @@ def _remove_constraints(self, constraints): super(Model, self)._remove_constraints(constraints) for constraint in constraints: self.problem.constraints.remove(constraint.name) - self.problem.clean() + self.problem.prune() def _get_variable_indices(self, names): vmap = dict(zip(self.variables, range(len(self.variables)))) diff --git a/src/optlang/osqp_interface.py b/src/optlang/osqp_interface.py index 3bfa5d0..b0c4c69 100644 --- a/src/optlang/osqp_interface.py +++ b/src/optlang/osqp_interface.py @@ -91,14 +91,21 @@ def osqp_settings(self): def solve(self): """Solve the OSQP problem.""" settings = self.osqp_settings() - P, q, A, bounds, _, _ = self.build(add_variable_constraints=True) + sp = self.build(add_variable_constraints=True) solver = osqp.OSQP() - if P is None: + if sp.P is None: # see https://github.com/cvxgrp/cvxpy/issues/898 settings.update({"adaptive_rho": 0, "rho": 1.0, "alpha": 1.0}) - solver.setup(P=P, q=q, A=A, l=bounds[:, 0], u=bounds[:, 1], **settings) # noqa + solver.setup( + P=sp.P, + q=sp.q, + A=sp.A, + l=sp.bounds[:, 0], + u=sp.bounds[:, 1], + **settings + ) if self._solution is not None: - if self.still_valid(A, bounds): + if self.still_valid(sp.A, sp.bounds): solver.warm_start(x=self._solution["x"], y=self._solution["y"]) if "rho" in self._solution: solver.update_settings(rho=self._solution["rho"]) @@ -109,7 +116,7 @@ def solve(self): self.primals = dict(zip(self.variables, solution.x)) self.vduals = dict(zip(self.variables, solution.y[nc : (nc + nv)])) if nc > 0: - self.cprimals = dict(zip(self.constraints, A.dot(solution.x)[0:nc])) + self.cprimals = dict(zip(self.constraints, sp.A.dot(solution.x)[0:nc])) self.duals = dict(zip(self.constraints, solution.y[0:nc])) if not np.isnan(solution.info.obj_val): self.obj_value = solution.info.obj_val * self.direction @@ -140,27 +147,22 @@ def still_valid(self, A, bounds): return valid -@six.add_metaclass(inheritdocstring) class Variable(mi.Variable): pass -@six.add_metaclass(inheritdocstring) class Constraint(mi.Constraint): _INDICATOR_CONSTRAINT_SUPPORT = False -@six.add_metaclass(inheritdocstring) class Objective(mi.Objective): pass -@six.add_metaclass(inheritdocstring) class Configuration(mi.Configuration): pass -@six.add_metaclass(inheritdocstring) class Model(mi.Model): ProblemClass = OSQPProblem status_map = _STATUS_MAP From 9c1f9706dd73e6f4012a73d1a43b2bc77cab359f Mon Sep 17 00:00:00 2001 From: Christian Diener Date: Wed, 13 Sep 2023 12:35:33 -0700 Subject: [PATCH 06/10] second round of reviews --- src/optlang/matrix_interface.py | 83 ++++++++++++++++++++++----------- 1 file changed, 55 insertions(+), 28 deletions(-) diff --git a/src/optlang/matrix_interface.py b/src/optlang/matrix_interface.py index e0bd6d7..b0c42a6 100644 --- a/src/optlang/matrix_interface.py +++ b/src/optlang/matrix_interface.py @@ -34,20 +34,19 @@ import abc import logging import pickle -from collections import defaultdict, namedtuple -from numpy import array, concatenate, Infinity, zeros - -from optlang import interface, symbolics, available_solvers -from optlang.expression_parsing import parse_optimization_expression -from optlang.exceptions import SolverError +from collections import defaultdict +from typing import NamedTuple +from numpy import Infinity, array, concatenate, zeros from scipy.sparse import csc_matrix +from optlang import available_solvers, interface, symbolics +from optlang.exceptions import SolverError +from optlang.expression_parsing import parse_optimization_expression from optlang.symbolics import add, mul log = logging.getLogger(__name__) - _STATUS_MAP = defaultdict(lambda: interface.UNDEFINED) _LP_METHODS = ("auto", ) @@ -56,11 +55,41 @@ _TYPES = ("continuous", "binary", "integer") -SparseProblem = namedtuple( - "SparseProblem", - ["P", "q", "A", "bounds", "vbounds", "integer"] -) -"""A representation of the Problem in standard form.""" + +class SparseProblem(NamedTuple): + """A representation of the convex optimizatgion problem in standard form. + + This defines the problem in the form. + + ..math:: + + \text{minimize }& \frac{1}{2}x^T\mathbf{P}x + q^Tx \\ + \text{s.t.: }& bounds_{.0} \leq \mathbf{A} \leq bounds_{.1} \\ + & vbounds_{.0} \leq x \leq vbounds_{.1} \\ + & \{x_k \in \mathbb{N}^0 \text{ if } integer_k = 1 \} + + Attributes + ---------- + P : csc_matrix + A semidefinite positive matrix specifying the coefficients for the quadratic + objective. + q : array + A vector specfiying the linear objective coefficients. + A : csc_matrix + The constraint matrix. + bounds : array + The lower and upper bounds corresponding to the rows in `A`. + vbounds : array + The lower and upper bounds of the variables `x`. + integer : array + Indicator vector denoting if `x[k]` is integer. + """ + P : csc_matrix + q : array + A : csc_matrix + bounds : array + vbounds : array + integer : array class MatrixProblem(abc.ABC): @@ -72,7 +101,6 @@ class MatrixProblem(abc.ABC): """ def __init__(self, **kwargs): - super().__init__(**kwargs) self.variables = set() self.integer_vars = set() self.constraints = set() @@ -342,7 +370,7 @@ def rename_variable(self, old, new): class Variable(interface.Variable): def __init__(self, name, *args, **kwargs): - super(Variable, self).__init__(name=name, **kwargs) + super().__init__(name=name, **kwargs) @interface.Variable.type.setter def type(self, value): @@ -376,8 +404,8 @@ def name(self, value): class Constraint(interface.Constraint): _INDICATOR_CONSTRAINT_SUPPORT = False - def __init__(self, expression, sloppy=False, *args, **kwargs): - super(Constraint, self).__init__(expression, *args, sloppy=sloppy, **kwargs) + def __init__(self, expression, sloppy=False, **kwargs): + super().__init__(expression=expression, sloppy=sloppy, **kwargs) def set_linear_coefficients(self, coefficients): if self.problem is not None: @@ -471,16 +499,16 @@ def __iadd__(self, other): if self.problem is not None: problem_reference = self.problem self.problem._remove_constraint(self) - super(Constraint, self).__iadd__(other) + super().__iadd__(other) problem_reference._add_constraint(self, sloppy=False) else: - super(Constraint, self).__iadd__(other) + super().__iadd__(other) return self class Objective(interface.Objective): def __init__(self, expression, sloppy=False, **kwargs): - super(Objective, self).__init__(expression, sloppy=sloppy, **kwargs) + super().__init__(expression=expression, sloppy=sloppy, **kwargs) self._expression_expired = False if not (sloppy or self.is_Linear or self.is_Quadratic): raise ValueError( @@ -561,10 +589,9 @@ def __init__( verbosity=0, timeout=None, qp_method="auto", - *args, **kwargs ): - super(Configuration, self).__init__(*args, **kwargs) + super().__init__(**kwargs) self.lp_method = lp_method self.presolve = presolve self.verbosity = verbosity @@ -708,7 +735,7 @@ def _initialize_model_from_problem(self, problem, vc_mapping=None, offset=0): ub=self.problem.variable_ubs[name], problem=self, ) - super(Model, self)._add_variables([var]) + super()._add_variables([var]) for name in self.problem.constraints: # Since constraint expressions are lazily retrieved from the @@ -722,7 +749,7 @@ def _initialize_model_from_problem(self, problem, vc_mapping=None, offset=0): problem=self, ) - super(Model, self)._add_constraints([constr], sloppy=True) + super()._add_constraints([constr], sloppy=True) if vc_mapping is None: for constr in self.constraints: @@ -847,7 +874,7 @@ def _set_variable_bounds_on_problem(self, var_lb, var_ub): self.problem.variable_ubs[var.name] = float(ub) def _add_variables(self, variables): - super(Model, self)._add_variables(variables) + super()._add_variables(variables) self.problem.reset() for variable in variables: lb = -Infinity if variable.lb is None else float(variable.lb) @@ -874,7 +901,7 @@ def _remove_variables(self, variables): self.problem.prune() def _add_constraints(self, constraints, sloppy=False): - super(Model, self)._add_constraints(constraints, sloppy=sloppy) + super()._add_constraints(constraints, sloppy=sloppy) self.problem.reset() for constraint in constraints: # This needs to be done in order to not trigger constraint._get_expression() @@ -905,7 +932,7 @@ def _add_constraints(self, constraints, sloppy=False): ) def _remove_constraints(self, constraints): - super(Model, self)._remove_constraints(constraints) + super()._remove_constraints(constraints) for constraint in constraints: self.problem.constraints.remove(constraint.name) self.problem.prune() @@ -955,10 +982,10 @@ def from_lp(self, lp_problem_str): mod = cplex_interface.Model.from_lp(lp_problem_str) mod.configuration.lp_method = "auto" mod.configuration.qp_method = "auto" - return super(Model, self).clone(mod) + return super().clone(mod) else: from optlang import glpk_interface mod = glpk_interface.Model.from_lp(lp_problem_str) mod.configuration.lp_method = "auto" - return super(Model, self).clone(mod) + return super().clone(mod) From 1524a3281921e5cea04ce4bd66fdb4d54b297b52 Mon Sep 17 00:00:00 2001 From: Christian Diener Date: Wed, 13 Sep 2023 15:39:48 -0700 Subject: [PATCH 07/10] fix variable types and move direction dual hacks into the MatrixProblem class --- src/optlang/matrix_interface.py | 49 +++++++++++++++++++++++---------- src/optlang/osqp_interface.py | 18 +++++------- 2 files changed, 42 insertions(+), 25 deletions(-) diff --git a/src/optlang/matrix_interface.py b/src/optlang/matrix_interface.py index b0c42a6..3a93c50 100644 --- a/src/optlang/matrix_interface.py +++ b/src/optlang/matrix_interface.py @@ -366,6 +366,9 @@ def rename_variable(self, old, new): } self.variable_lbs[new] = self.variable_lbs.pop(old) self.variable_ubs[new] = self.variable_ubs.pop(old) + if old in self.integer_vars: + self.integer_vars.remove(old) + self.integer_vars.add(new) class Variable(interface.Variable): @@ -381,6 +384,15 @@ def type(self, value): + "The following variable types are available: " + ", ".join(_TYPES) ) + if value == "binary": + self.lb = 0 + self.ub = 1 + self.problem.problem.integer_vars.add(self.name) + elif value == "integer": + self.problem.problem.integer_vars.add(self.name) + elif value == "continuous": + if self.name in self.problem.problem.integer_vars: + self.problem.problem.integer_vars.remove(self.name) super(Variable, Variable).type.fset(self, value) def _get_primal(self): @@ -388,6 +400,8 @@ def _get_primal(self): @property def dual(self): + if self.problem.is_integer: + raise ValueError("Dual values are not well-defined for integer problems") if self.problem is not None: return self.problem.problem.vduals.get(self.name, None) return None @@ -462,11 +476,11 @@ def primal(self): @property def dual(self): + if self.problem.is_integer: + raise ValueError("Dual values are not well-defined for integer problems") if self.problem is None: return None d = self.problem.problem.duals.get(self.name, None) - if d is not None: - d = -d return d @interface.Constraint.name.setter @@ -582,6 +596,9 @@ def get_linear_coefficients(self, variables): class Configuration(interface.MathematicalProgrammingConfiguration): + lp_methods = _LP_METHODS + qp_methods = _QP_METHODS + def __init__( self, lp_method="auto", @@ -605,14 +622,14 @@ def __init__( @property def lp_method(self): """The algorithm used to solve LP problems.""" - return "auto" + return self.problem.problem.settings["lp_method"] @lp_method.setter def lp_method(self, lp_method): - if lp_method not in _LP_METHODS: + if lp_method not in self.lp_methods: raise ValueError( "LP Method %s is not valid (choose one of: %s)" - % (lp_method, ", ".join(_LP_METHODS)) + % (lp_method, ", ".join(self.lp_methods)) ) def _set_presolve(self, value): @@ -679,14 +696,14 @@ def __setstate__(self, state): @property def qp_method(self): """Change the algorithm used to optimize QP problems.""" - return "auto" + return self.problem.problem.settings["qp_method"] @qp_method.setter def qp_method(self, value): if value not in _QP_METHODS: raise ValueError( "%s is not a valid qp_method. Choose between %s" - % (value, str(_QP_METHODS)) + % (value, str(self.qp_methods)) ) self._qp_method = value @@ -698,10 +715,10 @@ def _set_feasibility(self, value): self.problem.problem.settings["dual_inf_tolerance"] = value def _get_integrality(self): - return 1e-6 + return self.problem.problem.settings["mip_tolerance"] def _set_integrality(self, value): - pass + self.problem.problem.settings["mip_tolerance"] = value def _get_optimality(self): return self.problem.problem.settings["optimality_tolerance"] @@ -782,7 +799,7 @@ def _initialize_model_from_problem(self, problem, vc_mapping=None, offset=0): linear_expression + quadratic_expression + offset, problem=self, direction={-1: "max", 1: "min"}[self.problem.direction], - name="osqp_objective", + name="matrix_objective", ) @property @@ -835,21 +852,25 @@ def _get_primal_values(self): return primal_values def _get_reduced_costs(self): - if len(self.problem.duals) == 0: + if self.is_integer: + raise ValueError("Dual values are not well-defined for integer problems") + if len(self.problem.vduals) == 0: raise SolverError("The problem has not been solved yet!") reduced_costs = [self.problem.vduals[v.name] for v in self._variables] return reduced_costs def _get_constraint_values(self): - if len(self.problem.primals) == 0: + if len(self.problem.cprimals) == 0: raise SolverError("The problem has not been solved yet!") constraint_primals = [self.problem.cprimals[c.name] for c in self._constraints] return constraint_primals def _get_shadow_prices(self): - if len(self.problem.primals) == 0: + if self.is_integer: + raise ValueError("Dual values are not well-defined for integer problems") + if len(self.problem.duals) == 0: raise SolverError("The problem has not been solved yet!") - dual_values = [-self.problem.duals[c.name] for c in self._constraints] + dual_values = [self.problem.duals[c.name] for c in self._constraints] return dual_values @property diff --git a/src/optlang/osqp_interface.py b/src/optlang/osqp_interface.py index b0c4c69..7b9c1ca 100644 --- a/src/optlang/osqp_interface.py +++ b/src/optlang/osqp_interface.py @@ -20,12 +20,10 @@ Make sure that 'import osqp' runs without error. """ import logging -import six import numpy as np import optlang.matrix_interface as mi from optlang import interface -from optlang.util import inheritdocstring log = logging.getLogger(__name__) @@ -57,7 +55,7 @@ "maximum iterations reached": interface.ITERATION_LIMIT, "unsolved": interface.SPECIAL, "problem non convex": interface.SPECIAL, - "non-existing-status": "Here for testing that missing statuses are handled." + "non-existing-status": "Here for testing that missing statuses are handled.", } @@ -91,18 +89,14 @@ def osqp_settings(self): def solve(self): """Solve the OSQP problem.""" settings = self.osqp_settings() + d = float(self.direction) sp = self.build(add_variable_constraints=True) solver = osqp.OSQP() if sp.P is None: # see https://github.com/cvxgrp/cvxpy/issues/898 settings.update({"adaptive_rho": 0, "rho": 1.0, "alpha": 1.0}) solver.setup( - P=sp.P, - q=sp.q, - A=sp.A, - l=sp.bounds[:, 0], - u=sp.bounds[:, 1], - **settings + P=sp.P, q=sp.q, A=sp.A, l=sp.bounds[:, 0], u=sp.bounds[:, 1], **settings ) if self._solution is not None: if self.still_valid(sp.A, sp.bounds): @@ -116,8 +110,10 @@ def solve(self): self.primals = dict(zip(self.variables, solution.x)) self.vduals = dict(zip(self.variables, solution.y[nc : (nc + nv)])) if nc > 0: - self.cprimals = dict(zip(self.constraints, sp.A.dot(solution.x)[0:nc])) - self.duals = dict(zip(self.constraints, solution.y[0:nc])) + self.cprimals = dict( + zip(self.constraints, sp.A.dot(solution.x)[0:nc] * d) + ) + self.duals = dict(zip(self.constraints, solution.y[0:nc] * d)) if not np.isnan(solution.info.obj_val): self.obj_value = solution.info.obj_val * self.direction self.status = solution.info.status From 0dc3a12386004c4afd1080d0aaec7ee8e424d56a Mon Sep 17 00:00:00 2001 From: Christian Diener Date: Thu, 14 Sep 2023 11:38:10 -0700 Subject: [PATCH 08/10] upodate release notes --- CHANGELOG.rst | 3 +++ 1 file changed, 3 insertions(+) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index e387d87..d35cdb2 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -3,6 +3,9 @@ Next Release ----- +* refactor the OSQP interface into a generic matrix interface to allow easy addition + of new solvers that expect an immutable problem in standard form as input + 1.7.0 ----- * remove deprecated numpy type casts From 04acdc2162b6e85b6d87c146029dbb07e819b478 Mon Sep 17 00:00:00 2001 From: Christian Diener Date: Thu, 14 Sep 2023 14:13:45 -0700 Subject: [PATCH 09/10] fix python version req, enable 3.11 --- .github/workflows/main.yml | 2 +- setup.cfg | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index a8d9d90..7a3cda6 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -20,7 +20,7 @@ jobs: fail-fast: false matrix: os: [ubuntu-latest, macos-latest, windows-latest] - python-version: ["3.8", "3.9", "3.10"] + python-version: ["3.8", "3.9", "3.10", "3.11"] steps: - uses: actions/checkout@v2 diff --git a/setup.cfg b/setup.cfg index 08c5254..e9c66b8 100644 --- a/setup.cfg +++ b/setup.cfg @@ -37,7 +37,7 @@ install_requires = six >=1.9 swiglpk >=5.0.8 sympy >=1.12.0 -python_requires = >=2.7, !=3.0.*, !=3.1.*, !=3.2.*, !=3.3.*, !=3.4.* +python_requires = >=3.8 tests_require = tox packages = find: From ab256d1fa30975a54b334e5210ea1336755febd0 Mon Sep 17 00:00:00 2001 From: Christian Diener Date: Thu, 14 Sep 2023 14:16:34 -0700 Subject: [PATCH 10/10] really enable 3.11 --- tox.ini | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tox.ini b/tox.ini index fee43a4..6392954 100644 --- a/tox.ini +++ b/tox.ini @@ -1,5 +1,5 @@ [tox] -envlist = isort, black, flake8, safety, docs, py3{8,9,10}, py3{8,9,10}-symengine +envlist = isort, black, flake8, safety, docs, py3{8,9,10,11}, py3{8,9,10,11}-symengine [gh-actions] python =