From bfb71ac3b1bc94b815476d09ea54a2576eb1f41f Mon Sep 17 00:00:00 2001 From: Kristian Jensen Date: Mon, 21 Aug 2017 14:57:28 +0200 Subject: [PATCH] 1.2.2 (#129) - Fixed a couple of gurobi interface bugs - Fixed a bug where cplex took forever to set a large objective (e.g. 10000 vars) - Experimental support for symengine (requires Python >= 3.5, install with `conda install -c symengine python-symengine`) --- optlang/cplex_interface.py | 48 +++++---- optlang/duality.py | 11 +- optlang/expression_parsing.py | 5 +- optlang/glpk_interface.py | 29 ++--- optlang/gurobi_interface.py | 101 +++++++++++------- optlang/interface.py | 71 ++++--------- optlang/scipy_interface.py | 11 +- optlang/symbolics.py | 130 +++++++++++++++++++++++ optlang/tests/abstract_test_cases.py | 57 +++++++--- optlang/tests/test_cplex_interface.py | 26 +++-- optlang/tests/test_expression_parsing.py | 32 ++++-- optlang/tests/test_glpk_interface.py | 73 ++++++++++--- optlang/tests/test_gurobi_interface.py | 26 +++-- optlang/tests/test_interface.py | 20 +++- optlang/tests/test_scipy_interface.py | 18 +++- optlang/util.py | 24 +++-- 16 files changed, 484 insertions(+), 198 deletions(-) create mode 100644 optlang/symbolics.py diff --git a/optlang/cplex_interface.py b/optlang/cplex_interface.py index 01492cb8..88b7669a 100644 --- a/optlang/cplex_interface.py +++ b/optlang/cplex_interface.py @@ -29,22 +29,18 @@ import os from six.moves import StringIO -import sympy -from sympy.core.add import _unevaluated_Add -from sympy.core.mul import _unevaluated_Mul -from sympy.core.singleton import S import cplex from cplex.exceptions import CplexSolverError from optlang import interface +from optlang import symbolics from optlang.util import inheritdocstring, TemporaryFilename from optlang.expression_parsing import parse_optimization_expression from optlang.exceptions import SolverError log = logging.getLogger(__name__) -Zero = S.Zero -One = S.One +from optlang.symbolics import add, mul, One, Zero _STATUS_MAP = { 'MIP_abort_feasible': interface.ABORTED, @@ -242,8 +238,8 @@ def _get_expression(self): cplex_problem = self.problem.problem cplex_row = cplex_problem.linear_constraints.get_rows(self.name) variables = self.problem._variables - expression = sympy.Add._from_args( - [sympy.Mul._from_args((sympy.RealNumber(cplex_row.val[i]), variables[ind])) for i, ind in + expression = add( + [mul((symbolics.Real(cplex_row.val[i]), variables[ind])) for i, ind in enumerate(cplex_row.ind)]) self._expression = expression return self._expression @@ -364,7 +360,7 @@ def _get_expression(self): if self.problem is not None and self._expression_expired and len(self.problem._variables) > 0: cplex_problem = self.problem.problem coeffs = cplex_problem.objective.get_linear() - expression = sympy.Add._from_args([coeff * var for coeff, var in zip(coeffs, self.problem._variables) if coeff != 0.]) + expression = add([coeff * var for coeff, var in zip(coeffs, self.problem._variables) if coeff != 0.]) if cplex_problem.objective.get_num_quadratic_nonzeros() > 0: expression += self.problem._get_quadratic_expression(cplex_problem.objective.get_quadratic()) self._expression = expression + getattr(self.problem, "_objective_offset", 0) @@ -616,9 +612,9 @@ def __init__(self, problem=None, *args, **kwargs): # lhs = _unevaluated_Add(*[val * variables[i - 1] for i, val in zip(row.ind, row.val)]) lhs = 0 if isinstance(lhs, int): - lhs = sympy.Integer(lhs) + lhs = symbolics.Integer(lhs) elif isinstance(lhs, float): - lhs = sympy.RealNumber(lhs) + lhs = symbolics.Real(lhs) if sense == 'E': constr = Constraint(lhs, lb=rhs, ub=rhs, name=name, problem=self) elif sense == 'G': @@ -650,10 +646,9 @@ def __init__(self, problem=None, *args, **kwargs): if 'CPLEX Error 1219:' not in str(e): raise e else: - linear_expression = _unevaluated_Add( - *[_unevaluated_Mul(sympy.RealNumber(coeff), variables[index]) for index, coeff in - enumerate(self.problem.objective.get_linear()) if coeff != 0.]) - + linear_expression = add( + [mul(symbolics.Real(coeff), variables[index]) for index, coeff in enumerate(self.problem.objective.get_linear()) if coeff != 0.] + ) try: quadratic = self.problem.objective.get_quadratic() except IndexError: @@ -716,7 +711,9 @@ def objective(self, value): if self._objective is not None: # Reset previous objective variables = self.objective.variables if len(variables) > 0: - self.problem.objective.set_linear([(variable.name, 0.) for variable in variables]) + name_list = [var.name for var in variables] + index_dict = {n: i for n, i in zip(name_list, self._get_variable_indices(name_list))} + self.problem.objective.set_linear([(index_dict[variable.name], 0.) for variable in variables]) if self.problem.objective.get_num_quadratic_variables() > 0: self.problem.objective.set_quadratic([0. for _ in range(self.problem.variables.get_num())]) super(Model, self.__class__).objective.fset(self, value) @@ -726,7 +723,9 @@ def objective(self, value): # self.problem.objective.set_offset(float(offset)) # Not available prior to 12.6.2 self._objective_offset = offset if linear_coefficients: - self.problem.objective.set_linear([var.name, float(coef)] for var, coef in linear_coefficients.items()) + name_list = [var.name for var in linear_coefficients] + index_dict = {n: i for n, i in zip(name_list, self._get_variable_indices(name_list))} + self.problem.objective.set_linear([index_dict[var.name], float(coef)] for var, coef in linear_coefficients.items()) for key, coef in quadratic_coeffients.items(): if len(key) == 1: @@ -905,4 +904,17 @@ def _get_quadratic_expression(self, quadratic=None): terms.append(0.5 * val * self._variables[i_name] ** 2) else: pass # Only look at upper triangle - return _unevaluated_Add(*terms) + return add(terms) + + def _get_variable_indices(self, names): + # Cplex does not keep an index of variable names + # Getting indices thus takes roughly quadratic time + # If many indices are required an alternate and faster method is used, where the full name list must only + # be traversed once + if len(names) < 1000: + return self.problem.variables.get_indices(names) + else: + name_set = set(names) + all_names = self.problem.variables.get_names() + indices = {n: i for i, n in enumerate(all_names) if n in name_set} + return [indices[n] for n in names] diff --git a/optlang/duality.py b/optlang/duality.py index f7a47495..4aeae1d0 100644 --- a/optlang/duality.py +++ b/optlang/duality.py @@ -13,8 +13,7 @@ # See the License for the specific language governing permissions and # limitations under the License. -from sympy import Add, Mul - +import optlang # This function is very complex. Should maybe be refactored def convert_linear_problem_to_dual(model, sloppy=False, infinity=None, maintain_standard_form=True, prefix="dual_", dual_model=None): # NOQA @@ -84,6 +83,8 @@ def convert_linear_problem_to_dual(model, sloppy=False, infinity=None, maintain_ if constraint.lb != 0: dual_objective[const_var] = sign * constraint.lb for variable, coef in constraint.expression.as_coefficients_dict().items(): + if variable == 1: # pragma: no cover # For symengine + continue coefficients.setdefault(variable.name, {})[const_var] = sign * coef else: if constraint.lb is not None: @@ -105,6 +106,8 @@ def convert_linear_problem_to_dual(model, sloppy=False, infinity=None, maintain_ coefficients_dict = {constraint.expression.args[1]: constraint.expression.args[0]} for variable, coef in coefficients_dict.items(): + if variable == 1: # pragma: no cover # For symengine + continue if constraint.lb is not None: coefficients.setdefault(variable.name, {})[lb_var] = -sign * coef if constraint.ub is not None: @@ -131,7 +134,7 @@ def convert_linear_problem_to_dual(model, sloppy=False, infinity=None, maintain_ # Add dual constraints from primal objective primal_objective_dict = model.objective.expression.as_coefficients_dict() for variable in model.variables: - expr = Add(*((coef * dual_var) for dual_var, coef in coefficients[variable.name].items())) + expr = optlang.symbolics.add([(coef * dual_var) for dual_var, coef in coefficients[variable.name].items()]) obj_coef = primal_objective_dict[variable] if maximization: const = model.interface.Constraint(expr, lb=obj_coef, name=prefix + variable.name) @@ -140,7 +143,7 @@ def convert_linear_problem_to_dual(model, sloppy=False, infinity=None, maintain_ dual_model.add(const) # Make dual objective - expr = Add(*((coef * dual_var) for dual_var, coef in dual_objective.items() if coef != 0)) + expr = optlang.symbolics.add([(coef * dual_var) for dual_var, coef in dual_objective.items() if coef != 0]) if maximization: objective = model.interface.Objective(expr, direction="min") else: diff --git a/optlang/expression_parsing.py b/optlang/expression_parsing.py index 9a923ce7..b83c68f6 100644 --- a/optlang/expression_parsing.py +++ b/optlang/expression_parsing.py @@ -13,6 +13,8 @@ # See the License for the specific language governing permissions and # limitations under the License. +from optlang.symbolics import One + def parse_optimization_expression(obj, linear=True, quadratic=False, expression=None, **kwargs): """ @@ -79,9 +81,10 @@ def _parse_linear_expression(expression, expanded=False, **kwargs): coefficients = {} else: raise ValueError("Expression {} seems to be invalid".format(expression)) + for var in coefficients: if not (var.is_Symbol): - if var == 1: + if var == One: constant = var offset = float(coefficients[var]) elif expanded: diff --git a/optlang/glpk_interface.py b/optlang/glpk_interface.py index 26e5d898..6eb2562b 100644 --- a/optlang/glpk_interface.py +++ b/optlang/glpk_interface.py @@ -30,13 +30,11 @@ import os import six -import sympy -from sympy.core.add import _unevaluated_Add -from sympy.core.mul import _unevaluated_Mul from optlang.util import inheritdocstring, TemporaryFilename from optlang.expression_parsing import parse_optimization_expression from optlang import interface +from optlang import symbolics log = logging.getLogger(__name__) @@ -52,7 +50,7 @@ glp_get_mat_row, glp_get_row_ub, glp_get_row_type, glp_get_row_lb, glp_get_row_name, glp_get_obj_coef, \ glp_get_obj_dir, glp_scale_prob, GLP_SF_AUTO, glp_get_num_int, glp_get_num_bin, glp_mip_col_val, \ glp_mip_obj_val, glp_mip_status, GLP_ETMLIM, glp_adv_basis, glp_read_lp, glp_mip_row_val, \ - get_col_primals, get_col_duals, get_row_primals, get_row_duals + get_col_primals, get_col_duals, get_row_primals, get_row_duals, glp_delete_prob @@ -165,8 +163,8 @@ def _get_expression(self): nnz = glp_get_mat_row(self.problem.problem, self._index, ia, da) constraint_variables = [self.problem._variables[glp_get_col_name(self.problem.problem, ia[i])] for i in range(1, nnz + 1)] - expression = sympy.Add._from_args( - [sympy.Mul._from_args((sympy.RealNumber(da[i]), constraint_variables[i - 1])) for i in + expression = symbolics.add( + [symbolics.mul((symbolics.Real(da[i]), constraint_variables[i - 1])) for i in range(1, nnz + 1)]) self._expression = expression return self._expression @@ -308,9 +306,9 @@ def term_generator(): for index in range(1, glp_get_num_cols(self.problem.problem) + 1): coeff = glp_get_obj_coef(self.problem.problem, index) if coeff != 0.: - yield (sympy.RealNumber(coeff), variables[index - 1]) + yield (symbolics.Real(coeff), variables[index - 1]) - expression = sympy.Add._from_args([sympy.Mul._from_args(term) for term in term_generator()]) + expression = symbolics.add([symbolics.mul(term) for term in term_generator()]) self._expression = expression + getattr(self.problem, "_objective_offset", 0) self._expression_expired = False return self._expression @@ -531,9 +529,9 @@ def __init__(self, problem=None, *args, **kwargs): ) log.exception() if isinstance(lhs, int): - lhs = sympy.Integer(lhs) + lhs = symbolics.Integer(lhs) elif isinstance(lhs, float): - lhs = sympy.RealNumber(lhs) + lhs = symbolics.Real(lhs) constraint_id = glp_get_row_name(self.problem, j) for variable in constraint_variables: try: @@ -551,9 +549,10 @@ def __init__(self, problem=None, *args, **kwargs): for index in range(1, glp_get_num_cols(problem) + 1) ) self._objective = Objective( - _unevaluated_Add( - *[_unevaluated_Mul(sympy.RealNumber(term[0]), term[1]) for term in term_generator if - term[0] != 0.]), + symbolics.add( + [symbolics.mul((symbolics.Real(term[0]), term[1])) for term in term_generator if + term[0] != 0.] + ), problem=self, direction={GLP_MIN: 'min', GLP_MAX: 'max'}[glp_get_obj_dir(self.problem)]) glp_scale_prob(self.problem, GLP_SF_AUTO) @@ -580,6 +579,10 @@ def __setstate__(self, repr_dict): if repr_dict['glpk_status'] == 'optimal': self.optimize() # since the start is an optimal solution, nothing will happen here + # def __del__(self): # To make sure that the glpk problem is deleted when this is garbage collected + # Gotcha: When objects with a __del__ method are part of a referencing cycle, the entire cycle is never automatically garbage collected + # glp_delete_prob(self.problem) + @property def objective(self): return self._objective diff --git a/optlang/gurobi_interface.py b/optlang/gurobi_interface.py index 9d5ddff5..2191be1c 100644 --- a/optlang/gurobi_interface.py +++ b/optlang/gurobi_interface.py @@ -27,11 +27,11 @@ import os import six +from functools import partial from optlang import interface from optlang.util import inheritdocstring, TemporaryFilename from optlang.expression_parsing import parse_optimization_expression -import sympy -from sympy.core.add import _unevaluated_Add +from optlang import symbolics import gurobipy @@ -211,10 +211,9 @@ def _get_expression(self): if internal_var_name == self.name + '_aux': continue variable = self.problem._variables[internal_var_name] - coeff = sympy.RealNumber(row.getCoeff(i)) - terms.append(sympy.Mul._from_args((coeff, variable))) - sympy.Add._from_args(terms) - self._expression = sympy.Add._from_args(terms) + coeff = symbolics.Real(row.getCoeff(i)) + terms.append(symbolics.mul((coeff, variable))) + self._expression = symbolics.add(terms) return self._expression def __str__(self): @@ -238,7 +237,14 @@ def problem(self, value): @property def primal(self): if self.problem is not None: - return self._internal_constraint.Slack + aux_var = self.problem.problem.getVarByName(self._internal_constraint.getAttr('ConstrName') + '_aux') + if aux_var is not None: + aux_val = aux_var.X + else: + aux_val = 0 + return (self._internal_constraint.RHS - + self._internal_constraint.Slack + + aux_val) # return self._round_primal_to_bounds(primal_from_solver) # Test assertions fail # return primal_from_solver else: @@ -327,9 +333,10 @@ def __init__(self, expression, sloppy=False, *args, **kwargs): @property def value(self): if getattr(self, "problem", None) is not None: + gurobi_problem = self.problem.problem try: - return self.problem.problem.getAttr("ObjVal") + getattr(self.problem, "_objective_offset", 0) - except gurobipy.GurobiError: # TODO: let this check the actual error message + return gurobi_problem.getAttr("ObjVal") + getattr(self.problem, "_objective_offset", 0) + except (gurobipy.GurobiError, AttributeError): # TODO: let this check the actual error message return None else: return None @@ -361,9 +368,10 @@ def _get_expression(self): if self.problem is not None and self._expression_expired and len(self.problem._variables) > 0: grb_obj = self.problem.problem.getObjective() terms = [] + variables = self.problem._variables for i in range(grb_obj.size()): - terms.append(grb_obj.getCoeff(i) * self.problem.variables[grb_obj.getVar(i).getAttr('VarName')]) - expression = sympy.Add._from_args(terms) + terms.append(grb_obj.getCoeff(i) * variables[grb_obj.getVar(i).getAttr('VarName')]) + expression = symbolics.add(terms) # TODO implement quadratic objectives self._expression = expression + getattr(self.problem, "_objective_offset", 0) self._expression_expired = False @@ -387,7 +395,8 @@ def lp_method(self): def lp_method(self, value): if value not in _LP_METHODS: raise ValueError("Invalid LP method. Please choose among: " + str(list(_LP_METHODS))) - self.problem.problem.params.Method = _LP_METHODS[value] + if getattr(self, "problem", None) is not None: + self.problem.problem.params.Method = _LP_METHODS[value] @property def presolve(self): @@ -395,12 +404,13 @@ def presolve(self): @presolve.setter def presolve(self, value): - if value is True: - self.problem.problem.params.Presolve = 2 - elif value is False: - self.problem.problem.params.Presolve = 0 - else: - self.problem.problem.params.Presolve = -1 + if getattr(self, "problem", None) is not None: + if value is True: + self.problem.problem.params.Presolve = 2 + elif value is False: + self.problem.problem.params.Presolve = 0 + else: + self.problem.problem.params.Presolve = -1 self._presolve = value @property @@ -429,20 +439,37 @@ def timeout(self, value): self.problem.problem.params.TimeLimit = value self._timeout = value + def __getstate__(self): + return {'presolve': self.presolve, 'verbosity': self.verbosity, 'timeout': self.timeout} + + def __setstate__(self, state): + self.__init__() + for key, val in six.iteritems(state): + setattr(self, key, val) + + def _get_feasibility(self): + return getattr(self.problem.problem.params, "FeasibilityTol") + + def _set_feasibility(self, value): + return setattr(self.problem.problem.params, "FeasibilityTol", value) + + def _get_optimality(self): + return getattr(self.problem.problem.params, "OptimalityTol") + + def _set_optimality(self, value): + return setattr(self.problem.problem.params, "OptimalityTol", value) + + def _get_integrality(self): + return getattr(self.problem.problem.params, "IntFeasTol") + + def _set_integrality(self, value): + return setattr(self.problem.problem.params, "IntFeasTol", value) + def _tolerance_functions(self): return { - "feasibility": ( - lambda: self.problem.problem.params.FeasibilityTol, - lambda x: setattr(self.problem.problem.params, "FeasibilityTol", x) - ), - "optimality": ( - lambda: self.problem.problem.params.OptimalityTol, - lambda x: setattr(self.problem.problem.params, "OptimalityTol", x) - ), - "integrality": ( - lambda: self.problem.problem.params.IntFeasTol, - lambda x: setattr(self.problem.problem.params, "IntFeasTol", x) - ) + "feasibility": (self._get_feasibility, self._set_feasibility), + "optimality": (self._get_optimality, self._set_optimality), + "integrality": (self._get_integrality, self._set_integrality) } @@ -477,9 +504,10 @@ def __init__(self, problem=None, *args, **kwargs): name = gurobi_constraint.getAttr("ConstrName") rhs = gurobi_constraint.RHS row = self.problem.getRow(gurobi_constraint) - lhs = _unevaluated_Add( - *[sympy.RealNumber(row.getCoeff(i)) * self.variables[row.getVar(i).VarName] for i in - range(row.size())]) + lhs = symbolics.add( + [symbolics.Real(row.getCoeff(i)) * self.variables[row.getVar(i).VarName] for i in + range(row.size())] + ) if sense == '=': constraint = Constraint(lhs, name=name, lb=rhs, ub=rhs, problem=self) @@ -493,9 +521,10 @@ def __init__(self, problem=None, *args, **kwargs): super(Model, self)._add_constraints(constraints, sloppy=True) gurobi_objective = self.problem.getObjective() - linear_expression = _unevaluated_Add( - *[sympy.RealNumber(gurobi_objective.getCoeff(i)) * self.variables[gurobi_objective.getVar(i).VarName] - for i in range(gurobi_objective.size())]) + linear_expression = symbolics.add( + [symbolics.Real(gurobi_objective.getCoeff(i)) * self.variables[gurobi_objective.getVar(i).VarName] + for i in range(gurobi_objective.size())] + ) self._objective = Objective( linear_expression, diff --git a/optlang/interface.py b/optlang/interface.py index a02cf115..0f7bd00c 100644 --- a/optlang/interface.py +++ b/optlang/interface.py @@ -29,17 +29,14 @@ import sys import uuid import warnings +import sympy import six from optlang.exceptions import IndicatorConstraintsNotSupported -import sympy - -from sympy.core.assumptions import _assume_rules -from sympy.core.facts import FactKB -from sympy.core.expr import Expr from optlang.util import parse_expr, expr_to_json, is_numeric, SolverTolerances +from optlang import symbolics from .container import Container @@ -77,7 +74,7 @@ # noinspection PyShadowingBuiltins -class Variable(sympy.Symbol): +class Variable(symbolics.Symbol): """Optimization variables. Variable objects are used to represents each variable of the optimization problem. When the optimization is @@ -148,46 +145,22 @@ def clone(cls, variable, **kwargs): """ return cls(variable.name, lb=variable.lb, ub=variable.ub, type=variable.type, **kwargs) - # def __new__(cls, name, **assumptions): - # - # if assumptions.get('zero', False): - # return S.Zero - # is_commutative = fuzzy_bool(assumptions.get('commutative', True)) - # if is_commutative is None: - # raise ValueError( - # '''Symbol commutativity must be True or False.''') - # assumptions['commutative'] = is_commutative - # for key in assumptions.keys(): - # assumptions[key] = bool(assumptions[key]) - # return sympy.Symbol.__xnew__(cls, name, uuid=str(int(round(1e16 * random.random()))), - # **assumptions) # uuid.uuid1() - - def __new__(cls, name, **kwargs): - if not isinstance(name, six.string_types): - raise TypeError("name should be a string, not %s" % repr(type(name))) - - obj = Expr.__new__(cls) - - obj.name = name - obj._assumptions = FactKB(_assume_rules) - obj._assumptions._tell('commutative', True) - obj._assumptions._tell('uuid', uuid.uuid1()) - - return obj - def __init__(self, name, lb=None, ub=None, type="continuous", problem=None, *args, **kwargs): # Ensure that name is str and not binary of unicode - some solvers only support string type in Python 2. if six.PY2: name = str(name) + if len(name) < 1: + raise ValueError('Variable name must not be empty string') + for char in name: if char.isspace(): raise ValueError( 'Variable names cannot contain whitespace characters. "%s" contains whitespace character "%s".' % ( name, char)) self._name = name - sympy.Symbol.__init__(name, *args, **kwargs) + symbolics.Symbol.__init__(self, name, *args, **kwargs) self._lb = lb self._ub = ub if self._lb is None and type == 'binary': @@ -457,13 +430,13 @@ def expression(self): @property def variables(self): """Variables in constraint/objective's expression.""" - return self.expression.atoms(sympy.Symbol) + return self.expression.atoms(Variable) def _canonicalize(self, expression): if isinstance(expression, float): - return sympy.RealNumber(expression) + return symbolics.Real(expression) elif isinstance(expression, int): - return sympy.Integer(expression) + return symbolics.Integer(expression) else: # expression = expression.expand() This would be a good way to canonicalize, but is quite slow return expression @@ -511,7 +484,12 @@ def is_Quadratic(self): is_quad = True return is_quad else: - poly = self.expression.as_poly(*self.variables) + if isinstance(self.expression, sympy.Basic): + sympy_expression = self.expression + else: + sympy_expression = sympy.sympify(self.expression) + # TODO: Find a way to do this with symengine (Poly is not part of symengine, 23 March 2017) + poly = sympy_expression.as_poly(*sympy_expression.atoms(sympy.Symbol)) if poly is None: return False else: @@ -888,7 +866,10 @@ def __eq__(self, other): def _canonicalize(self, expression): """For example, changes x + y to 1.*x + 1.*y""" expression = super(Objective, self)._canonicalize(expression) - expression *= 1. + if isinstance(expression, sympy.Basic): + expression *= 1. + else: # pragma: no cover # symengine + expression = (1. * expression).expand() return expression @property @@ -1180,16 +1161,10 @@ def objective(self): @objective.setter def objective(self, value): self.update() - # try: # Is this except ever needed? for atom in sorted(value.expression.atoms(Variable), key=lambda v: v.name): if isinstance(atom, Variable) and (atom.problem is None or atom.problem != self): self._pending_modifications.add_var.append(atom) self.update() - # except AttributeError as e: - # if isinstance(value.expression, six.types.FunctionType) or isinstance(value.expression, float): - # pass - # else: - # raise AttributeError(e) if self._objective is not None: self._objective.problem = None self._objective = value @@ -1323,7 +1298,7 @@ def _get_shadow_prices(self): def is_integer(self): return any(var.type in ("integer", "binary") for var in self.variables) - def __str__(self): + def __str__(self): # pragma: no cover if hasattr(self, "to_lp"): return self.to_lp() self.update() @@ -1507,9 +1482,9 @@ def _remove_variables(self, variables): replacements = dict([(variable, 0) for variable in variables]) for constraint_id in constraint_ids: constraint = self._constraints[constraint_id] - constraint._expression = constraint.expression.xreplace(replacements) + constraint._expression = constraint._expression.xreplace(replacements) if self.objective is not None: - self.objective._expression = self.objective.expression.xreplace(replacements) + self.objective._expression = self.objective._expression.xreplace(replacements) def _remove_variable(self, variable): self._remove_variables([variable]) diff --git a/optlang/scipy_interface.py b/optlang/scipy_interface.py index 52459b5d..d4ad1603 100644 --- a/optlang/scipy_interface.py +++ b/optlang/scipy_interface.py @@ -34,7 +34,7 @@ from optlang.expression_parsing import parse_optimization_expression from scipy.optimize import linprog from optlang.util import inheritdocstring -import sympy +from optlang import symbolics SCIPY_STATUS = { 0: interface.OPTIMAL, @@ -417,7 +417,7 @@ def ub(self, value): def coefficient_dict(self, names=True): if self.expression.is_Add: coefficient_dict = {variable: coef for variable, coef in - self.expression.as_coefficients_dict().items()} + self.expression.as_coefficients_dict().items() if variable.is_Symbol} elif self.expression.is_Atom and self.expression.is_Symbol: coefficient_dict = {self.expression: 1} elif self.expression.is_Mul and len(self.expression.args) <= 2: @@ -437,7 +437,7 @@ def set_linear_coefficients(self, coefficients): self.lb, self.ub = None, None coefficients_dict = self.coefficient_dict(names=False) coefficients_dict.update(coefficients) - self._expression = sympy.Add(*(v * k for k, v in coefficients_dict.items())) + self._expression = symbolics.add(*(v * k for k, v in coefficients_dict.items())) self.lb = lb self.ub = ub else: @@ -465,7 +465,7 @@ def _get_expression(self): coefficients_dict = { self.problem._variables[name]: coef for name, coef in coefficients_dict.items() if name in self.problem._variables } - self._expression = sympy.Add(*(v * k for k, v in coefficients_dict.items())) + self.problem.problem.offset + self._expression = symbolics.add(*(v * k for k, v in coefficients_dict.items())) + self.problem.problem.offset return self._expression @property @@ -484,7 +484,7 @@ def direction(self, value): def coefficient_dict(self): if self.expression.is_Add: coefficient_dict = {variable: coef for variable, coef in - self.expression.as_coefficients_dict().items()} + self.expression.as_coefficients_dict().items() if variable.is_Symbol} elif self.expression.is_Atom and self.expression.is_Symbol: coefficient_dict = {self.expression: 1} elif self.expression.is_Mul and len(self.expression.args) <= 2 and self.expression.args[1].is_Symbol: @@ -633,7 +633,6 @@ def objective(self, value): value.problem = self - if __name__ == "__main__": x1 = Variable('x1', lb=0) x2 = Variable('x2', lb=0) diff --git a/optlang/symbolics.py b/optlang/symbolics.py new file mode 100644 index 00000000..db0d9870 --- /dev/null +++ b/optlang/symbolics.py @@ -0,0 +1,130 @@ +# Copyright 2017 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. + +""" +This module contains a common interface to symbolic operations in Sympy and Symengine respectively. +All symbolic operations in the optlang codebase should use these functions. +""" + +from __future__ import division + +import os +import warnings +import six +import uuid + +# Read environment variable +SYMENGINE_PREFERENCE = os.environ.get("OPTLANG_USE_SYMENGINE", "") + +if SYMENGINE_PREFERENCE.lower() in ("false", "no", "off"): + USE_SYMENGINE = False +else: + try: + import symengine + import symengine.sympy_compat + from symengine.sympy_compat import Symbol as symengine_Symbol + except ImportError as e: + if SYMENGINE_PREFERENCE.lower() in ("true", "yes", "on"): + warnings.warn("Symengine could not be imported: " + str(e)) + USE_SYMENGINE = False + else: + USE_SYMENGINE = True + + +if USE_SYMENGINE: # pragma: no cover + import operator + from six.moves import reduce + + Integer = symengine.Integer + Real = symengine.RealDouble + Zero = Integer(0) + One = Integer(1) + NegativeOne = Integer(-1) + sympify = symengine.sympy_compat.sympify + + Add = symengine.Add + Mul = symengine.Mul + Pow = symengine.sympy_compat.Pow + + class Symbol(symengine_Symbol): + def __new__(cls, name, *args, **kwargs): + return symengine_Symbol.__new__(cls, name) + + def __init__(self, name, *args, **kwargs): + super(Symbol, self).__init__(name) + self._name = name + + def __repr__(self): + return self._name + + def __str__(self): + return self._name + + def __getnewargs__(self): + return (self._name, {}) + + def add(*args): + if len(args) == 1: + args = args[0] + return sum(args) + + def mul(*args): + if len(args) == 1: + args = args[0] + return reduce(operator.mul, args, 1) + +else: # Use sympy + import sympy + from sympy.core.assumptions import _assume_rules + from sympy.core.facts import FactKB + from sympy.core.expr import Expr + + Integer = sympy.Integer + Real = sympy.RealNumber + Zero = Integer(0) + One = Integer(1) + NegativeOne = Integer(-1) + sympify = sympy.sympify + + Add = sympy.Add + Mul = sympy.Mul + Pow = sympy.Pow + + class Symbol(sympy.Symbol): + def __new__(cls, name, **kwargs): + if not isinstance(name, six.string_types): + raise TypeError("name should be a string, not %s" % repr(type(name))) + + obj = sympy.Symbol.__new__(cls, str(uuid.uuid1())) + + obj.name = name + obj._assumptions = FactKB(_assume_rules) + obj._assumptions._tell('commutative', True) + obj._assumptions._tell('uuid', uuid.uuid1()) + + return obj + + def __init__(self, *args, **kwargs): + super(Symbol, self).__init__() + + def add(*args): + if len(args) == 1: + args = args[0] + return sympy.Add._from_args(args) + + def mul(*args): + if len(args) == 1: + args = args[0] + return sympy.Mul._from_args(args) diff --git a/optlang/tests/abstract_test_cases.py b/optlang/tests/abstract_test_cases.py index f8dd66d9..f5ce5ea2 100644 --- a/optlang/tests/abstract_test_cases.py +++ b/optlang/tests/abstract_test_cases.py @@ -63,6 +63,7 @@ def test_change_name(self): self.model.update() self.var.name = "test_2" self.assertEqual(self.var.name, "test_2") + self.assertEqual(str(self.var), self.var.name) self.model.remove(self.var) self.model.update() @@ -338,7 +339,7 @@ def test_create_empty_model(self): model = self.interface.Model() self.assertEqual(len(model.constraints), 0) self.assertEqual(len(model.variables), 0) - self.assertEqual(model.objective.expression, 0) + self.assertEqual(model.objective.expression - 0, 0) @abc.abstractmethod def test_pickle_ability(self): @@ -346,12 +347,12 @@ def test_pickle_ability(self): def test_pickle_empty_model(self): model = self.interface.Model() - self.assertEquals(model.objective.expression, 0) + self.assertEquals(model.objective.expression - 0, 0) self.assertEquals(len(model.variables), 0) self.assertEquals(len(model.constraints), 0) pickle_string = pickle.dumps(model) from_pickle = pickle.loads(pickle_string) - self.assertEquals(from_pickle.objective.expression, 0) + self.assertEquals(from_pickle.objective.expression - 0, 0) self.assertEquals(len(from_pickle.variables), 0) self.assertEquals(len(from_pickle.constraints), 0) @@ -442,7 +443,7 @@ def test_add_constraints(self): self.model.update() self.model.add(constr2) self.model.update() - self.model.add(constr3) + self.model.add(constr3, sloppy=True) self.model.update() self.model.add([constr4, constr5, constr6]) self.model.update() @@ -481,14 +482,14 @@ def test_objective_get_linear_coefficients(self): coefs = self.model.objective.get_linear_coefficients(self.model.variables) self.assertEqual(len(coefs), len(self.model.variables)) expr = sum(c * v for v, c in coefs.items()) - self.assertEqual(expr - self.model.objective.expression, 0) + self.assertEqual((expr - self.model.objective.expression).expand() - 0, 0) def test_constraint_get_linear_coefficients(self): constraint = self.model.constraints[5] coefs = constraint.get_linear_coefficients(self.model.variables) self.assertEqual(len(coefs), len(self.model.variables)) expr = sum(c * v for v, c in coefs.items()) - self.assertEqual(expr - constraint.expression, 0) + self.assertEqual((expr - constraint.expression).expand() - 0, 0) @abc.abstractmethod def test_change_of_constraint_is_reflected_in_low_level_solver(self): @@ -514,9 +515,8 @@ def test_change_variable_type(self): def test_change_constraint_bounds(self): pass - @abc.abstractmethod def test_initial_objective(self): - pass + self.assertEqual(self.model.objective.expression, 1.0 * self.model.variables["R_Biomass_Ecoli_core_w_GAM"]) def test_optimize(self): self.model.optimize() @@ -536,13 +536,23 @@ def test_optimize_milp(self): def test_change_objective(self): v1, v2 = self.model.variables.values()[0:2] self.model.objective = self.interface.Objective(1. * v1 + 1. * v2) - self.assertEqual(str(self.model.objective), 'Maximize\n1.0*R_PGK + 1.0*R_Biomass_Ecoli_core_w_GAM') + self.assertEqual(self.model.objective.direction, "max") + self.assertEqual( + (self.model.objective.expression - + (1.0 * self.model.variables["R_PGK"] + 1.0 * self.model.variables["R_Biomass_Ecoli_core_w_GAM"])).expand(), + 0. + ) self.model.objective = self.interface.Objective(v1 + v2) - self.assertEqual(str(self.model.objective), 'Maximize\n1.0*R_PGK + 1.0*R_Biomass_Ecoli_core_w_GAM') + self.assertEqual(self.model.objective.direction, "max") + self.assertEqual( + (self.model.objective.expression - + (1.0 * self.model.variables["R_PGK"] + 1.0 * self.model.variables["R_Biomass_Ecoli_core_w_GAM"])).expand(), + 0. + ) def test_number_objective(self): self.model.objective = self.interface.Objective(0.) - self.assertEqual(self.model.objective.expression, 0) + self.assertEqual(self.model.objective.expression - 0, 0) self.assertEqual(self.model.objective.direction, "max") self.assertEqual(self.model.optimize(), "optimal") @@ -667,13 +677,13 @@ def test_objective_set_linear_coefficients(self): self.assertAlmostEqual(y.primal, 0) obj.set_linear_coefficients({y: 1}) - self.assertEqual(obj.expression - (x + y), 0) + self.assertEqual((obj.expression - (x + y)).expand(), 0) self.assertEqual(model.optimize(), optlang.interface.OPTIMAL) self.assertAlmostEqual(x.primal, 2) self.assertAlmostEqual(y.primal, 2) obj.set_linear_coefficients({x: 0}) - self.assertEqual(obj.expression - y, 0) + self.assertEqual((obj.expression - y).expand(), 0) self.assertEqual(model.optimize(), optlang.interface.OPTIMAL) self.assertAlmostEqual(x.primal, 0) self.assertAlmostEqual(y.primal, 3) @@ -691,12 +701,12 @@ def test_constraint_set_linear_coefficients(self): self.assertAlmostEqual(x.primal, x.ub) c1.set_linear_coefficients({x: 1}) - self.assertEqual(c1.expression - (x + y), 0) + self.assertEqual((c1.expression - (x + y)).expand() - 0, 0) self.assertEqual(model.optimize(), optlang.interface.OPTIMAL) self.assertAlmostEqual(x.primal, 1) c1.set_linear_coefficients({x: 2}) - self.assertEqual(c1.expression - (2 * x + y), 0) + self.assertEqual((c1.expression - (2 * x + y)).expand() - 0, 0) self.assertEqual(model.optimize(), optlang.interface.OPTIMAL) self.assertAlmostEqual(x.primal, 0.5) @@ -721,14 +731,16 @@ def test_objective_expression_includes_constant(self): objective = self.model.objective self.model.objective = self.interface.Objective(objective.expression + 3) self.model.update() - self.assertEqual((self.model.objective.expression - (objective.expression + 3)).expand(), 0) + self.assertEqual((self.model.objective.expression - (objective.expression + 3.)).expand(), 0.) def test_is_integer(self): model = self.model self.assertFalse(model.is_integer) + self.assertFalse(optlang.interface.Model.is_integer.fget(model)) model.variables[0].type = "integer" self.assertTrue(model.is_integer) + self.assertTrue(optlang.interface.Model.is_integer.fget(model)) model.variables[0].type = "continuous" model.variables[1].type = "binary" @@ -799,6 +811,19 @@ def test_integer_batch_duals(self): self.assertEqual(model.reduced_costs[x.name], 0) self.assertEqual(model.shadow_prices[c.name], 1) + def test_large_objective(self): + model = self.interface.Model() + model.add([self.interface.Variable(str(i), lb=1) for i in range(1100)]) + model.optimize() + + obj = self.interface.Objective( + optlang.symbolics.add([optlang.symbolics.mul((optlang.symbolics.One, v)) for v in model.variables]), + direction="min" + ) + model.objective = obj + model.optimize() + self.assertAlmostEqual(model.objective.value, len(model.variables)) + @six.add_metaclass(abc.ABCMeta) class AbstractConfigurationTestCase(unittest.TestCase): diff --git a/optlang/tests/test_cplex_interface.py b/optlang/tests/test_cplex_interface.py index ae4871fb..48463b54 100644 --- a/optlang/tests/test_cplex_interface.py +++ b/optlang/tests/test_cplex_interface.py @@ -110,7 +110,7 @@ def test_get_primal(self): self.assertEqual(self.constraint.primal, None) self.model.optimize() self.assertEqual(self.model.status, 'optimal') - self.assertEqual(self.model.objective.value, 0.8739215069684305) + self.assertAlmostEqual(self.model.objective.value, 0.8739215069684305) print([constraint.primal for constraint in self.model.constraints]) for i, j in zip([constraint.primal for constraint in self.model.constraints], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, @@ -214,7 +214,12 @@ def test_change_of_constraint_is_reflected_in_low_level_solver(self): y = Variable('y', lb=-181133.3, ub=12000.) constraint = Constraint(0.3 * x + 0.4 * y, lb=-100, name='test') self.model.add(constraint) - self.assertEqual(self.model.constraints['test'].__str__(), 'test: -100 <= 0.4*y + 0.3*x') + self.assertEqual( + (self.model.constraints['test'].expression - (0.4 * y + 0.3 * x)).expand() - 0, + 0 + ) + + # self.assertEqual(self.model.constraints['test'].__str__(), 'test: -100 <= 0.4*y + 0.3*x') self.assertEqual(self.model.problem.linear_constraints.get_coefficients([('test', 'x'), ('test', 'y')]), [0.3, 0.4]) z = Variable('z', lb=3, ub=4, type='integer') @@ -222,8 +227,10 @@ def test_change_of_constraint_is_reflected_in_low_level_solver(self): self.assertEqual( self.model.problem.linear_constraints.get_coefficients([('test', 'x'), ('test', 'y'), ('test', 'z')]), [0.3, 0.4, 77.]) - self.assertEqual(self.model.constraints['test'].__str__(), 'test: -100 <= 0.4*y + 0.3*x + 77.0*z') - print(self.model) + self.assertEqual( + (self.model.constraints['test'].expression - (0.4 * y + 0.3 * x + 77.0 * z)).expand() - 0, + 0 + ) def test_constraint_set_problem_to_None_caches_the_latest_expression_from_solver_instance(self): x = Variable('x', lb=-83.3, ub=1324422.) @@ -233,7 +240,10 @@ def test_constraint_set_problem_to_None_caches_the_latest_expression_from_solver z = Variable('z', lb=2, ub=5, type='integer') constraint += 77. * z self.model.remove(constraint) - self.assertEqual(constraint.__str__(), 'test: -100 <= 0.4*y + 0.3*x + 77.0*z') + self.assertEqual(constraint.name, "test") + self.assertEqual(constraint.lb, -100) + self.assertEqual(constraint.ub, None) + self.assertEqual((constraint.expression - (0.4 * y + 0.3 * x + 77.0 * z)).expand(), 0) def test_change_of_objective_is_reflected_in_low_level_solver(self): x = Variable('x', lb=-83.3, ub=1324422.) @@ -298,9 +308,6 @@ def test_change_constraint_bounds(self): self.assertEqual(self.model.problem.linear_constraints.get_senses(constraint.name), "E") self.assertEqual(self.model.problem.linear_constraints.get_range_values(constraint.name), 0) - def test_initial_objective(self): - self.assertEqual(self.model.objective.expression.__str__(), '1.0*R_Biomass_Ecoli_core_w_GAM') - def test_iadd_objective(self): v2, v3 = self.model.variables.values()[1:3] self.model.objective += 2. * v2 - 3. * v3 @@ -330,7 +337,8 @@ def test_imul_objective(self): def test_set_copied_objective(self): obj_copy = copy.copy(self.model.objective) self.model.objective = obj_copy - self.assertEqual(self.model.objective.__str__(), 'Maximize\n1.0*R_Biomass_Ecoli_core_w_GAM') + self.assertEqual(self.model.objective.direction, "max") + self.assertEqual(self.model.objective.expression, 1.0 * self.model.variables["R_Biomass_Ecoli_core_w_GAM"]) def test_timeout(self): self.model.configuration.timeout = 0 diff --git a/optlang/tests/test_expression_parsing.py b/optlang/tests/test_expression_parsing.py index 2a0f458b..2b307b82 100644 --- a/optlang/tests/test_expression_parsing.py +++ b/optlang/tests/test_expression_parsing.py @@ -7,6 +7,20 @@ import unittest +def _quad_terms_to_expression(terms): + args = [] + for term, coef in terms.items(): + term = list(term) + args.append(coef * term[0] * term[-1]) + return sum(args) + + +def _compare_term_dicts(test_case, dict1, dict2): + for term, coef1 in dict1.items(): + coef2 = dict2[term] + test_case.assertEqual(coef1 - coef2, 0) + + class ExpressionParsingTestCase(unittest.TestCase): def setUp(self): self.vars = [Variable(name) for name in ["x", "y", "z", "v", "w"]] @@ -22,8 +36,8 @@ def test_parse_linear_expression(self): self.assertEqual(offset_const, 0) self.assertEqual(offset_obj, offset) - self.assertEqual(linear_terms_const, target) - self.assertEqual(linear_terms_obj, target) + _compare_term_dicts(self, linear_terms_const, target) + _compare_term_dicts(self, linear_terms_obj, target) self.assertEqual(quad_terms_const, {}) self.assertEqual(quad_terms_obj, {}) @@ -41,8 +55,10 @@ def test_parse_quadratic_expression(self): self.assertEqual(offset_obj, offset) self.assertEqual(linear_terms_const, {}) self.assertEqual(linear_terms_obj, {}) - self.assertEqual(quad_terms_const, target) - self.assertEqual(quad_terms_obj, target) + _compare_term_dicts(self, quad_terms_const, target) + _compare_term_dicts(self, quad_terms_obj, target) + self.assertEqual((_quad_terms_to_expression(quad_terms_obj) - (expr - offset)).expand(), 0) + self.assertEqual((_quad_terms_to_expression(quad_terms_const) - (expr - offset)).expand(), 0) def test_parse_non_expanded_quadratic_expression(self): x, y, z = self.vars[:3] @@ -57,7 +73,7 @@ def test_parse_non_expanded_quadratic_expression(self): self.assertEqual(offset_const, -4) self.assertEqual(offset_obj, -4 + offset) - self.assertEqual(linear_terms_const, linear_target) - self.assertEqual(linear_terms_obj, linear_target) - self.assertEqual(quad_terms_const, target) - self.assertEqual(quad_terms_obj, target) + _compare_term_dicts(self, linear_terms_const, linear_target) + _compare_term_dicts(self, linear_terms_obj, linear_target) + _compare_term_dicts(self, quad_terms_const, target) + _compare_term_dicts(self, quad_terms_obj, target) diff --git a/optlang/tests/test_glpk_interface.py b/optlang/tests/test_glpk_interface.py index acb99435..0b00e18d 100644 --- a/optlang/tests/test_glpk_interface.py +++ b/optlang/tests/test_glpk_interface.py @@ -91,7 +91,7 @@ def test_get_primal(self): 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.5546369406238147e-17, 0.0, -5.080374405378186e-29, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]): - self.assertAlmostEqual(i, j) + self.assertAlmostEqual(i, j, places=6) class ObjectiveTestCase(abstract_test_cases.AbstractObjectiveTestCase): @@ -162,14 +162,33 @@ def test_add_non_cplex_conform_variable(self): self.assertEqual(var_from_pickle.name, glp_get_col_name(repickled.problem, var_from_pickle._index)) def test_glpk_remove_variable(self): - var = self.model.variables.values()[0] - self.assertEqual(self.model.constraints['M_atp_c'].__str__(), - 'M_atp_c: 0.0 <= -1.0*R_ACKr - 1.0*R_ADK1 + 1.0*R_ATPS4r - 1.0*R_PGK - 1.0*R_SUCOAS - 59.81*R_Biomass_Ecoli_core_w_GAM - 1.0*R_GLNS - 1.0*R_GLNabc - 1.0*R_PFK - 1.0*R_PPCK - 1.0*R_PPS + 1.0*R_PYK - 1.0*R_ATPM <= 0.0') # noqa: E501 + var = self.model.variables[0] + self.assertEqual( + (self.model.constraints["M_atp_c"].expression - ( + -1.0 * self.model.variables.R_ACKr - 1.0 * self.model.variables.R_ADK1 + 1.0 * self.model.variables.R_ATPS4r - + 1.0 * self.model.variables.R_PGK - 1.0 * self.model.variables.R_SUCOAS - 59.81 * self.model.variables.R_Biomass_Ecoli_core_w_GAM - + 1.0 * self.model.variables.R_GLNS - 1.0 * self.model.variables.R_GLNabc - 1.0 * self.model.variables.R_PFK - + 1.0 * self.model.variables.R_PPCK - 1.0 * self.model.variables.R_PPS + 1.0 * self.model.variables.R_PYK - + 1.0 * self.model.variables.R_ATPM + )).expand() - 0, 0 + ) + self.assertEqual(self.model.constraints["M_atp_c"].lb, 0) + self.assertEqual(self.model.constraints["M_atp_c"].ub, 0) + self.assertEqual(var.problem, self.model) self.model.remove(var) self.model.update() - self.assertEqual(self.model.constraints['M_atp_c'].__str__(), - 'M_atp_c: 0.0 <= -1.0*R_ACKr - 1.0*R_ADK1 + 1.0*R_ATPS4r - 1.0*R_PGK - 1.0*R_SUCOAS - 1.0*R_GLNS - 1.0*R_GLNabc - 1.0*R_PFK - 1.0*R_PPCK - 1.0*R_PPS + 1.0*R_PYK - 1.0*R_ATPM <= 0.0') # noqa: E501 + + self.assertEqual( + (self.model.constraints["M_atp_c"].expression - ( + -1.0 * self.model.variables.R_ACKr - 1.0 * self.model.variables.R_ADK1 + 1.0 * self.model.variables.R_ATPS4r - + 1.0 * self.model.variables.R_SUCOAS - 1.0 * self.model.variables.R_PGK - + 1.0 * self.model.variables.R_GLNS - 1.0 * self.model.variables.R_GLNabc - 1.0 * self.model.variables.R_PFK - + 1.0 * self.model.variables.R_PPCK - 1.0 * self.model.variables.R_PPS + 1.0 * self.model.variables.R_PYK - + 1.0 * self.model.variables.R_ATPM + )).expand() - 0, 0 + ) + self.assertNotIn(var, self.model.variables.values()) self.assertEqual(glp_find_col(self.model.problem, var.name), 0) self.assertEqual(var.problem, None) @@ -255,14 +274,22 @@ def test_change_of_constraint_is_reflected_in_low_level_solver(self): constraint = Constraint(0.3 * x + 0.4 * y, lb=-100, name='test') self.assertEqual(constraint._index, None) self.model.add(constraint) - self.assertEqual(self.model.constraints['test'].__str__(), 'test: -100 <= 0.4*y + 0.3*x') + self.assertEqual( + (self.model.constraints["test"].expression - (0.4 * y + 0.3 * x)).expand(), 0 + ) + self.assertEqual(self.model.constraints["test"].lb, -100) + self.assertEqual(constraint._index, 73) z = Variable('z', lb=3, ub=10, type='integer') self.assertEqual(z._index, None) constraint += 77. * z self.assertEqual(z._index, 98) - self.assertEqual(self.model.constraints['test'].__str__(), 'test: -100 <= 0.4*y + 0.3*x + 77.0*z') - print(self.model) + + self.assertEqual( + (self.model.constraints["test"].expression - (0.4 * y + 0.3 * x + 77.0 * z)).expand(), 0 + ) + self.assertEqual(self.model.constraints["test"].lb, -100) + self.assertEqual(constraint._index, 73) def test_constraint_set_problem_to_None_caches_the_latest_expression_from_solver_instance(self): @@ -273,14 +300,23 @@ def test_constraint_set_problem_to_None_caches_the_latest_expression_from_solver z = Variable('z', lb=2, ub=5, type='integer') constraint += 77. * z self.model.remove(constraint) - self.assertEqual(constraint.__str__(), 'test: -100 <= 0.4*y + 0.3*x + 77.0*z') + self.assertEqual( + (constraint.expression - (0.4 * y + 0.3 * x + 77.0 * z)).expand() - 0, 0 + ) + self.assertEqual(constraint.lb, -100) + self.assertEqual(constraint.ub, None) def test_change_of_objective_is_reflected_in_low_level_solver(self): x = Variable('x', lb=-83.3, ub=1324422.) y = Variable('y', lb=-181133.3, ub=12000.) objective = Objective(0.3 * x + 0.4 * y, name='test', direction='max') self.model.objective = objective - self.assertEqual(self.model.objective.__str__(), 'Maximize\n0.4*y + 0.3*x') + + self.assertEqual( + (self.model.objective.expression - (0.4 * y + 0.3 * x)).expand() - 0, 0 + ) + self.assertEqual(self.model.objective.direction, "max") + self.assertEqual(glp_get_obj_coef(self.model.problem, x._index), 0.3) self.assertEqual(glp_get_obj_coef(self.model.problem, y._index), 0.4) for i in range(1, glp_get_num_cols(self.model.problem) + 1): @@ -288,7 +324,12 @@ def test_change_of_objective_is_reflected_in_low_level_solver(self): self.assertEqual(glp_get_obj_coef(self.model.problem, i), 0) z = Variable('z', lb=4, ub=4, type='integer') self.model.objective += 77. * z - self.assertEqual(self.model.objective.__str__(), 'Maximize\n0.4*y + 0.3*x + 77.0*z') + + self.assertEqual( + (self.model.objective.expression - (0.4 * y + 0.3 * x + 77.0 * z)).expand() - 0, 0 + ) + self.assertEqual(self.model.objective.direction, "max") + self.assertEqual(glp_get_obj_coef(self.model.problem, x._index), 0.3) self.assertEqual(glp_get_obj_coef(self.model.problem, y._index), 0.4) self.assertEqual(glp_get_obj_coef(self.model.problem, z._index), 77.) @@ -328,9 +369,6 @@ def test_change_constraint_bounds(self): self.assertNotEqual(inner_problem_bounds, inner_problem_bounds_new) self.assertEqual(bounds_new, inner_problem_bounds_new) - def test_initial_objective(self): - self.assertEqual(self.model.objective.expression.__str__(), '1.0*R_Biomass_Ecoli_core_w_GAM') - def test_iadd_objective(self): v2, v3 = self.model.variables.values()[1:3] self.model.objective += 2. * v2 - 3. * v3 @@ -363,7 +401,10 @@ def test_imul_objective(self): def test_set_copied_objective(self): obj_copy = copy.copy(self.model.objective) self.model.objective = obj_copy - self.assertEqual(self.model.objective.__str__(), 'Maximize\n1.0*R_Biomass_Ecoli_core_w_GAM') + self.assertEqual( + (self.model.objective.expression - 1.0 * self.model.variables.R_Biomass_Ecoli_core_w_GAM).expand() - 0, 0 + ) + self.assertEqual(self.model.objective.direction, "max") def test_timeout(self): self.model.configuration.timeout = 0 diff --git a/optlang/tests/test_gurobi_interface.py b/optlang/tests/test_gurobi_interface.py index 56a68128..cd5cdc3a 100644 --- a/optlang/tests/test_gurobi_interface.py +++ b/optlang/tests/test_gurobi_interface.py @@ -286,12 +286,20 @@ def test_change_of_constraint_is_reflected_in_low_level_solver(self): constraint = Constraint(0.3 * x + 0.4 * y, lb=-100, name='test') self.assertEqual(constraint._internal_constraint, None) self.model.add(constraint) - self.assertEqual(self.model.constraints['test'].__str__(), 'test: -100 <= 0.4*y + 0.3*x') + self.assertEqual(self.model.constraints['test'].lb, -100) + self.assertEqual( + (self.model.constraints['test'].expression - (0.4 * y + 0.3 * x)).expand() - 0, + 0 + ) z = Variable('z', lb=3, ub=10, type='integer') self.assertEqual(z._internal_variable, None) constraint += 77. * z self.assertEqual(z._internal_variable, self.model.problem.getVarByName('z')) - self.assertEqual(self.model.constraints['test'].__str__(), 'test: -100 <= 0.4*y + 0.3*x + 77.0*z') + self.assertEqual(self.model.constraints['test'].lb, -100) + self.assertEqual( + (self.model.constraints['test'].expression - (0.4 * y + 0.3 * x + 77.0 * z)).expand() - 0, + 0 + ) def test_constraint_set_problem_to_None_caches_the_latest_expression_from_solver_instance(self): x = Variable('x', lb=-83.3, ub=1324422.) @@ -301,7 +309,10 @@ def test_constraint_set_problem_to_None_caches_the_latest_expression_from_solver z = Variable('z', lb=2, ub=5, type='integer') constraint += 77. * z self.model.remove(constraint) - self.assertEqual(constraint.__str__(), 'test: -100 <= 0.4*y + 0.3*x + 77.0*z') + self.assertEqual(constraint.lb, -100) + self.assertEqual( + (constraint.expression - (0.4 * y + 0.3 * x + 77.0 * z)).expand() - 0, 0 + ) def test_change_of_objective_is_reflected_in_low_level_solver(self): x = Variable('x', lb=-83.3, ub=1324422.) @@ -361,9 +372,6 @@ def test_change_constraint_bounds(self): print('bounds', bounds) self.assertEqual(bounds, inner_problem_bounds) - def test_initial_objective(self): - self.assertEqual(self.model.objective.expression.__str__(), '1.0*R_Biomass_Ecoli_core_w_GAM') - @unittest.skip('Not supported yet') def test_iadd_objective(self): v2, v3 = self.model.variables.values()[1:3] @@ -400,7 +408,11 @@ def test_imul_objective(self): def test_set_copied_objective(self): obj_copy = copy.copy(self.model.objective) self.model.objective = obj_copy - self.assertEqual(self.model.objective.__str__(), 'Maximize\n1.0*R_Biomass_Ecoli_core_w_GAM') + self.assertEqual(self.model.objective.direction, "max") + self.assertEqual( + (self.model.objective.expression - (1.0 * self.model.variables.R_Biomass_Ecoli_core_w_GAM)).expand() - 0, + 0 + ) def test_timeout(self): self.model.configuration.timeout = 0 diff --git a/optlang/tests/test_interface.py b/optlang/tests/test_interface.py index ed89f8a7..80da66a0 100644 --- a/optlang/tests/test_interface.py +++ b/optlang/tests/test_interface.py @@ -29,7 +29,7 @@ def test_create_empty_model(self): model = Model() self.assertEqual(len(model.constraints), 0) self.assertEqual(len(model.variables), 0) - self.assertEqual(model.objective.expression, 0) + self.assertEqual(model.objective.expression - 0, 0) def test_add_variable(self): var = Variable('z') @@ -100,11 +100,22 @@ def test_remove_variable_str(self): self.model.remove(var.name) self.assertNotIn(var, self.model.variables.values()) self.assertEqual(var.problem, None) - self.assertEqual(self.model.__str__(), 'Maximize\n2.0*x\nsubject to\nconstr1: 3 <= 1.0*x\nBounds\n0 <= x <= 10') + self.assertEqual(self.model.objective.direction, "max") + self.assertEqual( + (self.model.objective.expression - (2.0 * self.model.variables.x)).expand() - 0, 0 + ) + self.assertEqual(len(self.model.variables), 1) + self.assertEqual(len(self.model.constraints), 1) + self.assertEqual((self.model.variables[0].lb, self.model.variables[0].ub), (0, 10)) + self.assertEqual(self.model.constraints[0].lb, 3) + self.assertEqual( + (self.model.constraints[0].expression - (1.0 * self.model.variables.x)).expand() - 0, 0 + ) + # self.assertEqual(self.model.__str__(), 'Maximize\n2.0*x\nsubject to\nconstr1: 3 <= 1.0*x\nBounds\n0 <= x <= 10') def test_number_objective(self): self.model.objective = Objective(0.) - self.assertEqual(self.model.objective.__str__(), 'Maximize\n0') + self.assertIn('Maximize\n0', self.model.objective.__str__()) def test_add_differing_interface_type_raises(self): from optlang import glpk_interface as glpk @@ -167,6 +178,9 @@ class TestVariable(TestCase): def setUp(self): self.x = Variable("x") + def test_init_variable(self): + self.assertRaises(ValueError, Variable, '') + def test_set_wrong_bounds_on_binary_raises(self): self.assertRaises(ValueError, Variable, 'x', lb=-33, ub=0.3, type='binary') x = Variable('x', type='binary') diff --git a/optlang/tests/test_scipy_interface.py b/optlang/tests/test_scipy_interface.py index 48de39b2..ab600f54 100644 --- a/optlang/tests/test_scipy_interface.py +++ b/optlang/tests/test_scipy_interface.py @@ -200,7 +200,12 @@ def test_change_constraint_bounds(self): pass def test_initial_objective(self): - self.assertIn('1.0*BIOMASS_Ecoli_core_w_GAM', self.model.objective.expression.__str__(), ) + self.assertIn('BIOMASS_Ecoli_core_w_GAM', self.model.objective.expression.__str__(), ) + self.assertEqual( + (self.model.objective.expression - ( + 1.0 * self.model.variables.BIOMASS_Ecoli_core_w_GAM - + 1.0 * self.model.variables.BIOMASS_Ecoli_core_w_GAM_reverse_712e5)).expand() - 0, 0 + ) @unittest.skip("Not implemented yet") def test_iadd_objective(self): @@ -214,7 +219,10 @@ def test_imul_objective(self): def test_set_copied_objective(self): obj_copy = copy.copy(self.model.objective) self.model.objective = obj_copy - self.assertEqual(self.model.objective.__str__(), 'Maximize\n1.0*R_Biomass_Ecoli_core_w_GAM') + self.assertEqual(self.model.objective.direction, "max") + self.assertEqual( + (self.model.objective.expression - (1.0 * self.model.variables.R_Biomass_Ecoli_core_w_GAM)).expand() - 0, 0 + ) @unittest.skip("NA") def test_timeout(self): @@ -293,6 +301,9 @@ def test_scipy_coefficient_dict(self): model.objective = obj self.assertEqual(model.optimize(), optlang.interface.OPTIMAL) + def test_deepcopy(self): + self.skipTest("Not implemented yet") + def test_is_integer(self): self.skipTest("No integers with scipy") @@ -304,3 +315,6 @@ def test_integer_constraint_dual(self): def test_integer_batch_duals(self): self.skipTest("No duals with scipy") + + def test_large_objective(self): + self.skipTest("Quite slow and not necessary") diff --git a/optlang/util.py b/optlang/util.py index 4b83f11c..aafb3590 100644 --- a/optlang/util.py +++ b/optlang/util.py @@ -26,6 +26,8 @@ from subprocess import check_output from sympy.printing.str import StrPrinter import sympy +from optlang import symbolics +from optlang.symbolics import mul, add, Pow def solve_with_glpsol(glp_prob): @@ -194,19 +196,19 @@ def expr_to_json(expr): """ Converts a Sympy expression to a json-compatible tree-structure. """ - if isinstance(expr, sympy.Mul): + if isinstance(expr, symbolics.Mul): return {"type": "Mul", "args": [expr_to_json(arg) for arg in expr.args]} - elif isinstance(expr, sympy.Add): + elif isinstance(expr, symbolics.Add): return {"type": "Add", "args": [expr_to_json(arg) for arg in expr.args]} - elif isinstance(expr, sympy.Symbol): + elif isinstance(expr, symbolics.Symbol): return {"type": "Symbol", "name": expr.name} - elif isinstance(expr, sympy.Pow): + elif isinstance(expr, symbolics.Pow): return {"type": "Pow", "args": [expr_to_json(arg) for arg in expr.args]} elif isinstance(expr, (float, int)): return {"type": "Number", "value": expr} - elif isinstance(expr, sympy.Float): + elif isinstance(expr, symbolics.Real): return {"type": "Number", "value": float(expr)} - elif isinstance(expr, sympy.Integer): + elif isinstance(expr, symbolics.Integer): return {"type": "Number", "value": int(expr)} else: raise NotImplementedError("Type not implemented: " + str(type(expr))) @@ -222,18 +224,18 @@ def parse_expr(expr, local_dict=None): if local_dict is None: local_dict = {} if expr["type"] == "Add": - return sympy.Add._from_args([parse_expr(arg, local_dict) for arg in expr["args"]]) + return add([parse_expr(arg, local_dict) for arg in expr["args"]]) elif expr["type"] == "Mul": - return sympy.Mul._from_args([parse_expr(arg, local_dict) for arg in expr["args"]]) + return mul([parse_expr(arg, local_dict) for arg in expr["args"]]) elif expr["type"] == "Pow": - return sympy.Pow(parse_expr(arg, local_dict) for arg in expr["args"]) + return Pow(parse_expr(arg, local_dict) for arg in expr["args"]) elif expr["type"] == "Symbol": try: return local_dict[expr["name"]] except KeyError: - return sympy.Symbol(expr["name"]) + return symbolics.Symbol(expr["name"]) elif expr["type"] == "Number": - return sympy.sympify(expr["value"]) + return symbolics.sympify(expr["value"]) else: raise NotImplementedError(expr["type"] + " is not implemented")