From 5cfe2399e8579b3b5a3fc8a4d141a4ded09d1a6b Mon Sep 17 00:00:00 2001 From: Nikolaus Sonnenschein Date: Fri, 21 Oct 2016 12:31:52 +0200 Subject: [PATCH] Gurobi interface (#19) * feat: working gurobi interface * feat: add objective property * test: add test_change_of_objective_is_reflected_in_low_level_solver to gurobi_interface * fix: set proper verbosity at model initiation * test: add netlib tests for gurobi_interface * fix: get constr name * feat: implement set_linear_coefficients for gurobi_interface.Constraint * feat: implement set_linear_coefficients for gurobi_interface.Objective * feat: implement constraint.name property * feat: add _get_expression to gurobi_interface.Objective * fix: flake8 * fix: cplex -> gurobipy in conditional before TestMissingDependency --- optlang/cplex_interface.py | 2 +- optlang/gurobi_interface.py | 612 ++++++++++++++++++++++++++ optlang/interface.py | 4 +- tests/test_gurobi_solver.py | 572 ++++++++++++++++++++++++ tests/test_netlib_gurobi_interface.py | 168 +++++++ 5 files changed, 1356 insertions(+), 2 deletions(-) create mode 100644 optlang/gurobi_interface.py create mode 100644 tests/test_gurobi_solver.py create mode 100644 tests/test_netlib_gurobi_interface.py diff --git a/optlang/cplex_interface.py b/optlang/cplex_interface.py index a706b3c8..01ec0e51 100644 --- a/optlang/cplex_interface.py +++ b/optlang/cplex_interface.py @@ -339,7 +339,7 @@ def set_linear_coefficients(self, coefficients): self.problem.problem.objective.set_linear([(variable.name, float(coefficient)) for variable, coefficient in coefficients.items()]) self._expression_expired = True else: - raise Exception("Can't change coefficients if constraint is not associated with a model.") + raise Exception("Can't change coefficients if objective is not associated with a model.") @six.add_metaclass(inheritdocstring) diff --git a/optlang/gurobi_interface.py b/optlang/gurobi_interface.py new file mode 100644 index 00000000..586da387 --- /dev/null +++ b/optlang/gurobi_interface.py @@ -0,0 +1,612 @@ +# 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. + + +from warnings import warn + +warn("Be careful! The GUROBI interface is still under construction ...") + +import logging + +log = logging.getLogger(__name__) + +import six +import tempfile +from optlang import interface +from optlang.util import inheritdocstring +import sympy +from sympy.core.add import _unevaluated_Add + +import gurobipy + +_GUROBI_STATUS_TO_STATUS = { + gurobipy.GRB.LOADED: interface.LOADED, + gurobipy.GRB.OPTIMAL: interface.OPTIMAL, + gurobipy.GRB.INFEASIBLE: interface.INFEASIBLE, + gurobipy.GRB.INF_OR_UNBD: interface.INFEASIBLE_OR_UNBOUNDED, + gurobipy.GRB.UNBOUNDED: interface.UNBOUNDED, + gurobipy.GRB.CUTOFF: interface.CUTOFF, + gurobipy.GRB.ITERATION_LIMIT: interface.ITERATION_LIMIT, + gurobipy.GRB.NODE_LIMIT: interface.NODE_LIMIT, + gurobipy.GRB.TIME_LIMIT: interface.TIME_LIMIT, + gurobipy.GRB.SOLUTION_LIMIT: interface.SOLUTION_LIMIT, + gurobipy.GRB.INTERRUPTED: interface.INTERRUPTED, + gurobipy.GRB.NUMERIC: interface.NUMERIC, + gurobipy.GRB.SUBOPTIMAL: interface.SUBOPTIMAL, + gurobipy.GRB.INPROGRESS: interface.INPROGRESS +} + +_VTYPE_TO_GUROBI_VTYPE = {'continuous': gurobipy.GRB.CONTINUOUS, 'integer': gurobipy.GRB.INTEGER, + 'binary': gurobipy.GRB.BINARY} +_GUROBI_VTYPE_TO_VTYPE = dict((val, key) for key, val in _VTYPE_TO_GUROBI_VTYPE.items()) + + +def _constraint_lb_and_ub_to_gurobi_sense_rhs_and_range_value(lb, ub): + """Helper function used by Constraint and Model""" + if lb is None and ub is None: + raise Exception("Free constraint ...") + elif lb is None: + sense = '<' + rhs = float(ub) + range_value = 0. + elif ub is None: + sense = '>' + rhs = float(lb) + range_value = 0. + elif lb == ub: + sense = '=' + rhs = float(lb) + range_value = 0. + elif lb > ub: + raise ValueError("Lower bound is larger than upper bound.") + else: + sense = '=' + rhs = float(lb) + range_value = float(ub - lb) + return sense, rhs, range_value + + +@six.add_metaclass(inheritdocstring) +class Variable(interface.Variable): + def __init__(self, name, *args, **kwargs): + super(Variable, self).__init__(name, **kwargs) + self._original_name = name # due to bug with getVarByName ... TODO: file bug report on gurobi google group + + @property + def _internal_variable(self): + if self.problem: + + internal_variable = self.problem.problem.getVarByName(self._original_name) + + assert internal_variable is not None + # if internal_variable is None: + # raise Exception('Cannot find variable {}') + return internal_variable + else: + return None + + @interface.Variable.lb.setter + def lb(self, value): + super(Variable, self.__class__).lb.fset(self, value) + if self.problem: + return self._internal_variable.setAttr('LB', value) + + @interface.Variable.ub.setter + def ub(self, value): + super(Variable, self.__class__).ub.fset(self, value) + if self.problem: + return self._internal_variable.setAttr('UB', value) + + @interface.Variable.type.setter + def type(self, value): + if self.problem: + return self._internal_variable.setAttr('VType', _VTYPE_TO_GUROBI_VTYPE[value]) + super(Variable, Variable).type.fset(self, value) + + @property + def primal(self): + if self.problem: + return self._internal_variable.getAttr('X') + else: + return None + + @property + def dual(self): + if self.problem: + return self._internal_variable.getAttr('RC') + else: + return None + + @interface.Variable.name.setter + def name(self, value): + if getattr(self, 'problem', None) is not None: + self._internal_variable.setAttr('VarName', value) + self.problem.problem.update() + super(Variable, Variable).name.fset(self, value) + + +@six.add_metaclass(inheritdocstring) +class Constraint(interface.Constraint): + _INDICATOR_CONSTRAINT_SUPPORT = False + + def __init__(self, expression, *args, **kwargs): + super(Constraint, self).__init__(expression, *args, **kwargs) + if self.ub is not None and self.lb is not None and self.lb > self.ub: + raise ValueError( + "Lower bound %f is larger than upper bound %f in constraint %s" % + (self.lb, self.ub, self) + ) + + def set_linear_coefficients(self, coefficients): + if self.problem is not None: + grb_constraint = self.problem.problem.getConstrByName(self.name) + for var, coeff in six.iteritems(coefficients): + self.problem.problem.chgCoeff(grb_constraint, self.problem.problem.getVarByName(var.name), float(coeff)) + else: + raise Exception("Can't change coefficients if constraint is not associated with a model.") + + @property + def _internal_constraint(self): + if self.problem is not None: + return self.problem.problem.getConstrByName(self.name) + + def _get_expression(self): + if self.problem is not None: + gurobi_problem = self.problem.problem + gurobi_constraint = self._internal_constraint + row = gurobi_problem.getRow(gurobi_constraint) + terms = [] + for i in range(row.size()): + internal_var_name = row.getVar(i).VarName + 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) + return self._expression + + def __str__(self): + if self.problem is not None: + self.problem.problem.update() + return super(Constraint, self).__str__() + + @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 = None + else: + self._problem = value + + @property + def primal(self): + if self.problem is not None: + return self._internal_constraint.Slack + # return self._round_primal_to_bounds(primal_from_solver) # Test assertions fail + # return primal_from_solver + else: + return None + + @property + def dual(self): + if self.problem is not None: + return self._internal_constraint.Pi + else: + return None + + @interface.Constraint.name.setter + def name(self, value): + if getattr(self, 'problem', None) is not None: + grb_constraint = self.problem.problem.getConstrByName(self.name) + grb_constraint.setAttr('ConstrName', value) + self._name = value + + @interface.Constraint.lb.setter + def lb(self, value): + if getattr(self, 'problem', None) is not None: + if self.ub is not None and value > self.ub: + raise ValueError( + "Lower bound %f is larger than upper bound %f in constraint %s" % + (value, self.ub, self) + ) + gurobi_constraint = self.problem.problem.getConstrByName(self.name) + sense, rhs, range_value = _constraint_lb_and_ub_to_gurobi_sense_rhs_and_range_value(value, self.ub) + if range_value == 0.: + gurobi_constraint.setAttr('Sense', sense) + gurobi_constraint.setAttr('RHS', rhs) + else: + aux_var = self.problem.problem.getVarByName(gurobi_constraint.getAttr('ConstrName') + '_aux') + if aux_var is None: + aux_var = self.problem.problem.addVar(name=gurobi_constraint.getAttr('ConstrName') + '_aux', lb=0, ub=range_value) + row = self.problem.problem.getRow(gurobi_constraint) + updated_row = row - aux_var + self.problem.problem.remove(gurobi_constraint) + self.problem.problem.update() + self.problem.problem.addConstr(updated_row, sense, rhs, self.name) + else: + aux_var.setAttr("UB", range_value) + self._lb = value + + @interface.Constraint.ub.setter + def ub(self, value): + if getattr(self, 'problem', None) is not None: + if self.lb is not None and value < self.lb: + raise ValueError( + "Upper bound %f is less than lower bound %f in constraint %s" % + (value, self.lb, self) + ) + gurobi_constraint = self.problem.problem.getConstrByName(self.name) + sense, rhs, range_value = _constraint_lb_and_ub_to_gurobi_sense_rhs_and_range_value(self.lb, value) + if range_value == 0.: + gurobi_constraint.setAttr('Sense', sense) + gurobi_constraint.setAttr('RHS', rhs) + else: + aux_var = self.problem.problem.getVarByName(gurobi_constraint.getAttr('ConstrName') + '_aux') + if aux_var is None: + aux_var = self.problem.problem.addVar(name=gurobi_constraint.getAttr('ConstrName') + '_aux', lb=0, ub=range_value) + row = self.problem.problem.getRow(gurobi_constraint) + updated_row = row - aux_var + self.problem.problem.remove(gurobi_constraint) + self.problem.problem.update() + self.problem.problem.addConstr(updated_row, sense, rhs, self.name) + else: + aux_var.setAttr("UB", range_value) + self._ub = value + + def __iadd__(self, 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, *args, **kwargs): + super(Objective, self).__init__(*args, **kwargs) + self._expression_expired = False + + @property + def value(self): + try: + return self.problem.problem.getAttr("ObjVal") + except gurobipy.GurobiError: # TODO: let this check the actual error message + return None + + @interface.Objective.direction.setter + def direction(self, value): + super(Objective, Objective).direction.__set__(self, value) + if getattr(self, 'problem', None) is not None: + self.problem.problem.setAttr('ModelSense', {'min': 1, 'max': -1}[value]) + + def set_linear_coefficients(self, coefficients): + if self.problem is not None: + grb_obj = self.problem.problem.getObjective() + for var, coeff in six.iteritems(coefficients): + grb_var = self.problem.problem.getVarByName(var.name) + grb_obj.remove(grb_var) + grb_obj.addTerms(float(coeff), grb_var) + self._expression_expired = True + else: + raise Exception("Can't change coefficients if objective is not associated with a model.") + + 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 = [] + 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) + # TODO implement quadratic objectives + self._expression = expression + self._expression_expired = False + return self._expression + + +@six.add_metaclass(inheritdocstring) +class Configuration(interface.MathematicalProgrammingConfiguration): + def __init__(self, lp_method='primal', tolerance=1e-9, presolve=False, verbosity=0, timeout=None, + solution_target="auto", qp_method="primal", *args, **kwargs): + super(Configuration, self).__init__(*args, **kwargs) + self.lp_method = lp_method + self.tolerance = tolerance + self.presolve = presolve + self.verbosity = verbosity + self.timeout = timeout + self.solution_target = solution_target + self.qp_method = qp_method + + @property + def presolve(self): + return self._presolve + + @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 + self._presolve = value + + @property + def verbosity(self): + return self._verbosity + + @verbosity.setter + def verbosity(self, value): + if self.problem is not None: + if value == 0: + self.problem.problem.params.OutputFlag = 0 + elif 0 < value <= 3: + self.problem.problem.params.OutputFlag = 1 + self._verbosity = value + + @property + def timeout(self): + return self._timeout + + @timeout.setter + def timeout(self, value): + if self.problem is not None: + if value is None: + self.problem.problem.params.TimeLimit = gurobipy.GRB.INFINITY + else: + self.problem.problem.params.TimeLimit = value + self._timeout = value + + +class Model(interface.Model): + """docstring for Model""" + + def __init__(self, problem=None, *args, **kwargs): + + super(Model, self).__init__(*args, **kwargs) + + if problem is None: + self.problem = gurobipy.Model() + self.problem.params.OutputFlag = 0 + if self.name is not None: + self.problem.setAttr('ModelName', self.name) + self.problem.update() + + elif isinstance(problem, gurobipy.Model): + self.problem = problem + variables = [] + for gurobi_variable in self.problem.getVars(): + variables.append(Variable( + gurobi_variable.getAttr("VarName"), + lb=gurobi_variable.getAttr("lb"), + ub=gurobi_variable.getAttr("ub"), + problem=self, + type=_GUROBI_VTYPE_TO_VTYPE[gurobi_variable.getAttr("vType")] + )) + super(Model, self)._add_variables(variables) + constraints = [] + for gurobi_constraint in self.problem.getConstrs(): + sense = gurobi_constraint.Sense + 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())]) + + if sense == '=': + constraint = Constraint(lhs, name=name, lb=rhs, ub=rhs, problem=self) + elif sense == '>': + constraint = Constraint(lhs, name=name, lb=rhs, ub=None, problem=self) + elif sense == '<': + constraint = Constraint(lhs, name=name, lb=None, ub=rhs, problem=self) + else: + raise ValueError('{} is not a valid sense'.format(sense)) + constraints.append(constraint) + 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())]) + + self._objective = Objective( + linear_expression, + problem=self, + direction={1: 'min', -1: 'max'}[self.problem.getAttr('ModelSense')] + ) + else: + raise TypeError("Provided problem is not a valid Gurobi model.") + + self.configuration = Configuration(problem=self, verbosity=0) + + def __getstate__(self): + tmp_file = tempfile.mktemp(suffix=".lp") + self.problem.write(tmp_file) + lp = open(tmp_file, 'rb').read() + repr_dict = {'lp': lp, 'status': self.status, 'config': self.configuration} + return repr_dict + + def __setstate__(self, repr_dict): + tmp_file = tempfile.mktemp(suffix=".lp") + open(tmp_file, 'wb').write(repr_dict['lp']) + problem = gurobipy.read(tmp_file) + # if repr_dict['status'] == 'optimal': # TODO: uncomment this + # # turn off logging completely, get's configured later + # problem.set_error_stream(None) + # problem.set_warning_stream(None) + # problem.set_log_stream(None) + # problem.set_results_stream(None) + # problem.solve() # since the start is an optimal solution, nothing will happen here + self.__init__(problem=problem) + self.configuration = Configuration.clone(repr_dict['config'], problem=self) # TODO: make configuration work + + def __str__(self): + self.problem.update() + tmp_file = tempfile.mktemp(suffix=".lp") + self.problem.update() + self.problem.write(tmp_file) + cplex_form = open(tmp_file).read() + return cplex_form + + @property + def objective(self): + return self._objective + + @objective.setter + def objective(self, value): + super(Model, self.__class__).objective.fset(self, value) + expression = self._objective._expression + if isinstance(expression, float) or isinstance(expression, int) or expression.is_Number: + pass + else: + if expression.is_Mul: + terms = (expression,) + elif expression.is_Add: + terms = expression.args + else: + raise ValueError( + "Provided objective %s doesn't seem to be appropriate." % + self._objective) + + grb_terms = [] + for term in terms: + factors = term.args + coeff = factors[0] + vars = factors[1:] + assert len(vars) <= 2 + var1 = self.problem.getVarByName(vars[0].name) + if len(vars) == 2: + var2 = self.problem.getVarByName(vars[1].name) + grb_terms.append(coeff * var1 * var2) + else: + if vars[0].is_Symbol: + grb_terms.append(coeff * var1) + elif vars[0].is_Pow: + grb_terms.append(coeff * var1**2) + + + grb_expression = gurobipy.quicksum(grb_terms) + + self.problem.setObjective(grb_expression, {'min': gurobipy.GRB.MINIMIZE, 'max': gurobipy.GRB.MAXIMIZE}[value.direction]) + value.problem = self + + def update(self): + super(Model, self).update(callback=self.problem.update) + + def _optimize(self): + self.problem.optimize() + self._status = _GUROBI_STATUS_TO_STATUS[self.problem.getAttr("Status")] + return self.status + + def _add_variables(self, variables): + super(Model, self)._add_variables(variables) + for variable in variables: + if variable.lb is None: + lb = -gurobipy.GRB.INFINITY + else: + lb = variable.lb + if variable.ub is None: + ub = gurobipy.GRB.INFINITY + else: + ub = variable.ub + gurobi_var = self.problem.addVar( + name=variable.name, + lb=lb, + ub=ub, + vtype=_VTYPE_TO_GUROBI_VTYPE[variable.type] + ) + variable._internal_var = gurobi_var + + def _remove_variables(self, variables): + # Not calling parent method to avoid expensive variable removal from sympy expressions + self.objective._expression = self.objective.expression.xreplace({var: 0 for var in variables}) + for variable in variables: + internal_variable = variable._internal_variable + del self._variables_to_constraints_mapping[variable.name] + variable.problem = None + del self._variables[variable.name] + self.problem.remove(internal_variable) + + def _add_constraints(self, constraints, sloppy=False): + super(Model, self)._add_constraints(constraints, sloppy=sloppy) + for constraint in constraints: + if constraint.lb is None and constraint.ub is None: + raise ValueError("optlang does not support free constraints in the gurobi interface.") + self.problem.update() + constraint._problem = None + if constraint.is_Linear: + if constraint.expression.is_Add: + coeff_dict = constraint.expression.as_coefficients_dict() + variables = [self.problem.getVarByName(variable.name) for variable in coeff_dict.keys()] + coefficients = [float(val) for val in coeff_dict.values()] + elif constraint.expression.is_Mul: + variable = list(constraint.expression.atoms(sympy.Symbol))[0] + variables = [self.problem.getVarByName(variable.name)] + coefficients = [float(constraint.expression.coeff(variable))] + elif constraint.expression.is_Atom and constraint.expression.is_Symbol: + variables = [self.problem.getVarByName(constraint.expression.name)] + coefficients = [1.] + elif constraint.expression.is_Number: + variables = [] + coefficients = [] + else: + raise ValueError('Something is fishy with constraint %s' % constraint) + lhs = gurobipy.quicksum([coeff * var for coeff, var in zip(coefficients, variables)]) + sense, rhs, range_value = _constraint_lb_and_ub_to_gurobi_sense_rhs_and_range_value(constraint.lb, + constraint.ub) + if range_value != 0.: + aux_var = self.problem.addVar(name=constraint.name + '_aux', lb=0, ub=range_value) + self.problem.update() + lhs = lhs - aux_var + rhs = constraint.ub + self.problem.addConstr(lhs, sense, rhs, name=constraint.name) + else: + raise ValueError( + "GUROBI currently only supports linear constraints. %s is not linear." % self) + # self.problem.addQConstr(lhs, sense, rhs) + constraint.problem = self + + def _remove_constraints(self, constraints): + self.problem.update() + internal_constraints = [constraint._internal_constraint for constraint in constraints] + super(Model, self)._remove_constraints(constraints) + for internal_constraint in internal_constraints: + self.problem.remove(internal_constraint) + + +if __name__ == '__main__': + x = Variable('x', lb=0, ub=10) + y = Variable('y', lb=0, ub=10) + z = Variable('z', lb=-100, ub=99.) + constr = Constraint(0.3 * x + 0.4 * y + 66. * z, lb=-100, ub=0., name='test') + + solver = Model(problem=gurobipy.read("tests/data/model.lp")) + + solver.add(z) + solver.add(constr) + print(solver) + print(solver.optimize()) + print(solver.status) diff --git a/optlang/interface.py b/optlang/interface.py index 8ba01a28..76250c93 100644 --- a/optlang/interface.py +++ b/optlang/interface.py @@ -1065,13 +1065,14 @@ def remove(self, stuff): raise TypeError( "Cannot remove %s. It neither a variable or constraint." % stuff) - def update(self): + def update(self, callback=int): """Process all pending model modifications.""" # print(self._pending_modifications) add_var = self._pending_modifications.add_var if len(add_var) > 0: self._add_variables(add_var) self._pending_modifications.add_var = [] + callback() add_constr = self._pending_modifications.add_constr if len(add_constr) > 0: @@ -1094,6 +1095,7 @@ def update(self): if len(rm_var) > 0: self._remove_variables(rm_var) self._pending_modifications.rm_var = [] + callback() rm_constr = self._pending_modifications.rm_constr if len(rm_constr) > 0: diff --git a/tests/test_gurobi_solver.py b/tests/test_gurobi_solver.py new file mode 100644 index 00000000..badedd2f --- /dev/null +++ b/tests/test_gurobi_solver.py @@ -0,0 +1,572 @@ +# Copyright (c) 2016 Novo Nordisk Foundation Center for Biosustainability, DTU. +# See LICENSE for details. + +import unittest + +try: + import gurobipy +except ImportError as e: + + class TestMissingDependency(unittest.TestCase): + + @unittest.skip('Missing dependency - ' + str(e)) + def test_fail(self): + pass + +else: + + import copy + + import random + + import os + + import nose + import pickle + + from optlang.gurobi_interface import Variable, Constraint, Model, Objective + + random.seed(666) + TESTMODELPATH = os.path.join(os.path.dirname(__file__), 'data/model.lp') + TESTMILPMODELPATH = os.path.join(os.path.dirname(__file__), 'data/simple_milp.lp') + + + class VariableTestCase(unittest.TestCase): + def setUp(self): + self.var = Variable('test') + + def test_internal_variable(self): + self.assertEqual(self.var._internal_variable, None) + + def test_set_wrong_type_raises(self): + self.assertRaises(Exception, setattr, self.var, 'type', 'ketchup') + + def test_get_primal(self): + self.assertEqual(self.var.primal, None) + model = Model(problem=gurobipy.read(TESTMODELPATH)) + model.optimize() + for i, j in zip([var.primal for var in model.variables], [0.8739215069684306, -16.023526143167608, 16.023526143167604, -14.71613956874283, 14.71613956874283, 4.959984944574658, 4.959984944574657, 4.959984944574658, 3.1162689467973905e-29, 2.926716099010601e-29, 0.0, 0.0, -6.112235045340358e-30, -5.6659435396316186e-30, 0.0, -4.922925402711085e-29, 0.0, 9.282532599166613, 0.0, 6.00724957535033, 6.007249575350331, 6.00724957535033, -5.064375661482091, 1.7581774441067828, 0.0, 7.477381962160285, 0.0, 0.22346172933182767, 45.514009774517454, 8.39, 0.0, 6.007249575350331, 0.0, -4.541857463865631, 0.0, 5.064375661482091, 0.0, 0.0, 2.504309470368734, 0.0, 0.0, -22.809833310204958, 22.809833310204958, 7.477381962160285, 7.477381962160285, 1.1814980932459636, 1.496983757261567, -0.0, 0.0, 4.860861146496815, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 5.064375661482091, 0.0, 5.064375661482091, 0.0, 0.0, 1.496983757261567, 10.000000000000002, -10.0, 0.0, 0.0, 0.0, 0.0, 0.0, -29.175827135565804, 43.598985311997524, 29.175827135565804, 0.0, 0.0, 0.0, -1.2332237321082153e-29, 3.2148950476847613, 38.53460965051542, 5.064375661482091, 0.0, -1.2812714099825612e-29, -1.1331887079263237e-29, 17.530865429786694, 0.0, 0.0, 0.0, 4.765319193197458, -4.765319193197457, 21.79949265599876, -21.79949265599876, -3.2148950476847613, 0.0, -2.281503094067127, 2.6784818505075303, 0.0]): + self.assertAlmostEqual(i, j) + + def test_get_dual(self): + self.assertEqual(self.var.dual, None) + model = Model(problem=gurobipy.read(TESTMODELPATH)) + model.optimize() + print([var.dual for var in model.variables]) + for i, j in zip([var.dual for var in model.variables], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.022916186593776214, 0.0, 0.0, 0.0, -0.03437427989066433, 0.0, -0.007638728864592076, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.005092485909728051, 0.0, 0.0, 0.0, 0.0, -0.005092485909728049, 0.0, 0.0, -0.005092485909728071, 0.0, 0.0, 0.0, -0.06110983091673658, -0.005092485909728054, 0.0, -0.003819364432296038, -0.005092485909728044, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.03946676580039238, 0.0, 0.0, -0.005092485909728051, 0.0, -0.0012731214774320122, 0.0, -0.09166474637510488, 0.0, 0.0, 0.0, -0.0, -0.045832373187552435, 0.0, 0.0, -0.09166474637510488, -0.005092485909728051, -0.07002168125876065, 0.0, -0.06874855978132864, -0.0012731214774320126, 0.0, 0.0, 0.0, -0.0012731214774320161, -0.003819364432296038, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.04073988727782439, -0.045832373187552435, -0.0012731214774320083, 0.0, 0.0, 0.0, 0.0, 0.0, -0.03437427989066433, 0.0, 0.0, -0.04837861614241646]): + self.assertAlmostEqual(i, j) + + def test_setting_lower_bound_higher_than_upper_bound_raises(self): + model = Model(problem=gurobipy.read(TESTMODELPATH)) + self.assertRaises(ValueError, setattr, model.variables[0], 'lb', 10000000000.) + + def test_setting_nonnumerical_bounds_raises(self): + model = Model(problem=gurobipy.read(TESTMODELPATH)) + self.assertRaises(Exception, setattr, model.variables[0], 'lb', 'Chicken soup') + + def test_changing_variable_names_is_reflected_in_the_solver(self): + model = Model(problem=gurobipy.read(TESTMODELPATH)) + for i, variable in enumerate(model.variables): + print(variable._internal_variable is not None) + print(variable.problem.name) + variable.name = "var"+str(i) + print(variable.problem.name) + print(variable.name) + print(variable._internal_variable is not None) + self.assertEqual(variable.name, "var"+str(i)) + self.assertEqual(variable._internal_variable.getAttr('VarName'), "var"+str(i)) + + def test_setting_bounds(self): + model = Model(problem=gurobipy.read(TESTMODELPATH)) + var = model.variables[0] + var.lb = 1 + self.assertEqual(var.lb, 1) + model.problem.update() + self.assertEqual(var._internal_variable.getAttr('LB'), 1) + var.ub = 2 + self.assertEqual(var.ub, 2) + model.problem.update() + self.assertEqual(var._internal_variable.getAttr('UB'), 2) + + + class ConstraintTestCase(unittest.TestCase): + def setUp(self): + self.model = Model(problem=gurobipy.read(TESTMODELPATH)) + self.constraint = Constraint(Variable('chip') + Variable('chap'), name='woodchips', lb=100) + + def test_get_primal(self): + self.assertEqual(self.constraint.primal, None) + self.model.optimize() + 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, 4.048900234729145e-15, 0.0, 0.0, 0.0, -3.55971196577979e-16, 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, 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) + + def test_get_dual(self): + self.assertEqual(self.constraint.dual, None) + self.model.optimize() + for i, j in zip([constraint.dual for constraint in self.model.constraints], [-0.04710549466498445, -0.0420130087552564, -0.0420130087552564, -0.09166474637510486, -0.09039162489767284, -0.024189308071208226, -0.022916186593776214, -0.03437427989066433, -0.03437427989066433, -0.028008672503504264, -0.07129480273619268, -0.029281793980936277, 0.0, -0.06238295239416859, -0.06110983091673658, 0.005092485909728049, -0.005092485909728049, -0.07129480273619268, -0.0, -0.0, 0.0, -0.0521979805747125, -0.06747543830389663, -0.04073988727782439, -0.03946676580039238, -0.09803035376226493, -0.104395961149425, -0.0, -0.0, -0.09166474637510488, -0.04837861614241646, -0.045832373187552435, -0.0521979805747125, -0.09803035376226493, -0.09166474637510488, -0.0751141671684887, -0.07002168125876065, -0.07002168125876065, -0.06874855978132864, -0.019096822161480183, -0.0, -0.0, 0.0012731214774320122, -0.0, -0.07129480273619268, -0.042013008755256404, -0.04073988727782439, -0.04837861614241646, -0.045832373187552435, 0.007638728864592073, 0.0, 0.008911850342024082, 0.0, -0.0, -0.0, -0.0, -0.0, -0.0420130087552564, -0.0420130087552564, -0.0012731214774320122, -0.0, -0.03564740136809634, -0.03437427989066433, 0.0, -0.0025462429548640244, -0.08275289603308078, -0.08275289603308078, -0.11330781149144906, -0.050924859097280485, -0.04837861614241646, -0.054744223529576516, -0.08275289603308078]): + self.assertAlmostEqual(i, j) + + def test_change_constraint_name(self): + constraint = copy.copy(self.constraint) + self.assertEqual(constraint.name, 'woodchips') + constraint.name = 'ketchup' + self.assertEqual(constraint.name, 'ketchup') + self.assertEqual([constraint.name for constraint in self.model.constraints], ['M_13dpg_c', 'M_2pg_c', 'M_3pg_c', 'M_6pgc_c', 'M_6pgl_c', 'M_ac_c', 'M_ac_e', 'M_acald_c', 'M_acald_e', 'M_accoa_c', 'M_acon_C_c', 'M_actp_c', 'M_adp_c', 'M_akg_c', 'M_akg_e', 'M_amp_c', 'M_atp_c', 'M_cit_c', 'M_co2_c', 'M_co2_e', 'M_coa_c', 'M_dhap_c', 'M_e4p_c', 'M_etoh_c', 'M_etoh_e', 'M_f6p_c', 'M_fdp_c', 'M_for_c', 'M_for_e', 'M_fru_e', 'M_fum_c', 'M_fum_e', 'M_g3p_c', 'M_g6p_c', 'M_glc_D_e', 'M_gln_L_c', 'M_gln_L_e', 'M_glu_L_c', 'M_glu_L_e', 'M_glx_c', 'M_h2o_c', 'M_h2o_e', 'M_h_c', 'M_h_e', 'M_icit_c', 'M_lac_D_c', 'M_lac_D_e', 'M_mal_L_c', 'M_mal_L_e', 'M_nad_c', 'M_nadh_c', 'M_nadp_c', 'M_nadph_c', 'M_nh4_c', 'M_nh4_e', 'M_o2_c', 'M_o2_e', 'M_oaa_c', 'M_pep_c', 'M_pi_c', 'M_pi_e', 'M_pyr_c', 'M_pyr_e', 'M_q8_c', 'M_q8h2_c', 'M_r5p_c', 'M_ru5p_D_c', 'M_s7p_c', 'M_succ_c', 'M_succ_e', 'M_succoa_c', 'M_xu5p_D_c'] + ) + for i, constraint in enumerate(self.model.constraints): + constraint.name = 'c'+ str(i) + self.assertEqual([constraint.name for constraint in self.model.constraints], ['c' + str(i) for i in range(0, len(self.model.constraints))]) + + def test_setting_lower_bound_higher_than_upper_bound_raises(self): + model = Model(problem=gurobipy.read(TESTMODELPATH)) + self.assertRaises(ValueError, setattr, model.constraints[0], 'lb', 10000000000.) + # + # def test_setting_nonnumerical_bounds_raises(self): + # model = Model(problem=glpk_read_cplex(TESTMODELPATH)) + # self.assertRaises(Exception, setattr, model.constraints[0], 'lb', 'Chicken soup') + + def test_set_linear_coefficients(self): + constraint = self.model.constraints.M_atp_c + constraint.set_linear_coefficients({self.model.variables.R_Biomass_Ecoli_core_w_GAM: 666.}) + self.model.update() + row = self.model.problem.getRow(self.model.problem.getConstrByName(constraint.name)) + for i in range(row.size()): + col_name = row.getVar(i).getAttr('VarName') + if col_name == 'R_Biomass_Ecoli_core_w_GAM': + self.assertEqual(row.getCoeff(i), 666.) + + + class ObjectiveTestCase(unittest.TestCase): + def setUp(self): + problem = gurobipy.read(TESTMODELPATH) + self.model = Model(problem=problem) + self.obj = self.model.objective + + def test_change_direction(self): + self.obj.direction = "min" + self.assertEqual(self.obj.direction, "min") + self.assertEqual(self.model.problem.getAttr('ModelSense'), gurobipy.GRB.MAXIMIZE) + self.model.update() + self.assertEqual(self.model.problem.getAttr('ModelSense'), gurobipy.GRB.MINIMIZE) + + self.obj.direction = "max" + self.assertEqual(self.obj.direction, "max") + self.assertEqual(self.model.problem.getAttr('ModelSense'), gurobipy.GRB.MINIMIZE) + self.model.update() + self.assertEqual(self.model.problem.getAttr('ModelSense'), gurobipy.GRB.MAXIMIZE) + + def test_set_linear_objective_coefficients(self): + self.model.objective.set_linear_coefficients({self.model.variables.R_TPI: 666.}) + self.model.update() + grb_obj = self.model.problem.getObjective() + for i in range(grb_obj.size()): + if 'R_TPI' == grb_obj.getVar(i).getAttr('VarName'): + self.assertEqual(grb_obj.getCoeff(i), 666.) + + + class SolverTestCase(unittest.TestCase): + def setUp(self): + # problem = gurobipy.Model() + problem = gurobipy.read(TESTMODELPATH) + self.model = Model(problem=problem) + + def test_create_empty_model(self): + model = Model() + self.assertEqual(model.problem.getAttr('NumVars'), 0) + self.assertEqual(model.problem.getAttr('NumConstrs'), 0) + self.assertEqual(model.name, None) + self.assertEqual(model.problem.getAttr('ModelName'), '') + model = Model(name="empty_problem") + self.assertEqual(model.problem.getAttr('ModelName'), 'empty_problem') + + def test_pickle_ability(self): + self.model.optimize() + value = self.model.objective.value + pickle_string = pickle.dumps(self.model) + from_pickle = pickle.loads(pickle_string) + from_pickle.optimize() + self.assertAlmostEqual(value, from_pickle.objective.value) + self.assertEqual([(var.lb, var.ub, var.name, var.type) for var in from_pickle.variables.values()], + [(var.lb, var.ub, var.name, var.type) for var in self.model.variables.values()]) + self.assertEqual([(constr.lb, constr.ub, constr.name) for constr in from_pickle.constraints], + [(constr.lb, constr.ub, constr.name) for constr in self.model.constraints]) + + def test_copy(self): + self.model.optimize() + value = self.model.objective.value + model_copy = copy.copy(self.model) + self.assertNotEqual(id(self.model), id(model_copy)) + self.assertNotEqual(id(self.model.problem), id(model_copy.problem)) + model_copy.optimize() + self.assertAlmostEqual(value, model_copy.objective.value) + self.assertEqual([(var.lb, var.ub, var.name, var.type) for var in model_copy.variables.values()], + [(var.lb, var.ub, var.name, var.type) for var in self.model.variables.values()]) + self.assertEqual([(constr.lb, constr.ub, constr.name) for constr in model_copy.constraints], + [(constr.lb, constr.ub, constr.name) for constr in self.model.constraints]) + + def test_deepcopy(self): + model_copy = copy.deepcopy(self.model) + self.assertNotEqual(id(self.model), id(model_copy)) + self.assertNotEqual(id(self.model.problem), id(model_copy.problem)) + + def test_config_gets_copied_too(self): + self.assertEquals(self.model.configuration.verbosity, 0) + self.model.configuration.verbosity = 3 + model_copy = copy.copy(self.model) + self.assertEquals(model_copy.configuration.verbosity, 3) + + def test_init_from_existing_problem(self): + self.assertEqual(len(self.model.variables), len(self.model.problem.getVars())) + self.assertEqual(len(self.model.constraints), len(self.model.problem.getConstrs())) + self.assertEqual(self.model.variables.keys(), + [var.VarName for var in self.model.problem.getVars()]) + + self.assertEqual(self.model.constraints.keys(), + [constr.ConstrName for constr in self.model.problem.getConstrs()]) + + def test_add_variable(self): + var = Variable('x') + self.model.add(var) + print(self.model._pending_modifications) + self.assertTrue(var in self.model.variables.values()) + self.assertEqual(self.model.variables.values().count(var), 1) + self.assertEqual(self.model.variables['x'].problem, var.problem) + print(var.name) + print(self.model.problem.getVars()) + print(self.model._pending_modifications) + self.model.update() + print(self.model._pending_modifications) + print(self.model.problem.getVars()) + self.assertEqual(self.model.problem.getVarByName(var.name).getAttr('VType'), gurobipy.GRB.CONTINUOUS) + var = Variable('y', lb=-13) + self.model.add(var) + self.assertTrue(var in self.model.variables.values()) + self.model.problem.update() + self.assertEqual(self.model.problem.getVarByName(var.name).getAttr('VType'), gurobipy.GRB.CONTINUOUS) + self.assertEqual(self.model.variables['x'].lb, None) + self.assertEqual(self.model.variables['x'].ub, None) + self.assertEqual(self.model.variables['y'].lb, -13) + self.assertEqual(self.model.variables['x'].ub, None) + var = Variable('x_with_ridiculously_long_variable_name_asdffffffffasdfasdfasdfasdfasdfasdfasdf') + self.model.add(var) + self.assertTrue(var in self.model.variables.values()) + self.assertEqual(self.model.variables.values().count(var), 1) + # var = Variable('x_with_ridiculously_long_variable_name_asdffffffffasdfasdfasdfasdfasdfasdfasdf') + # self.assertRaises(Exception, self.model.add, var) + # self.assertEqual(len(self.model.variables), len(self.model.problem.getVars())) + + def test_add_integer_var(self): + var = Variable('int_var', lb=-13, ub=500, type='integer') + self.model.add(var) + self.assertEqual(self.model.variables['int_var'].type, 'integer') + self.assertEqual(self.model.problem.getVarByName(var.name).getAttr('VType'), gurobipy.GRB.INTEGER) + self.assertEqual(self.model.variables['int_var'].ub, 500) + self.assertEqual(self.model.variables['int_var'].lb, -13) + + def test_add_non_cplex_conform_variable(self): + var = Variable('12x!!@#5_3', lb=-666, ub=666) + self.model.add(var) + self.assertTrue(var in self.model.variables.values()) + self.model.problem.update() + self.assertEqual(var.name, self.model.problem.getVarByName(var.name).VarName) + self.assertEqual(self.model.variables['12x!!@#5_3'].lb, -666) + self.assertEqual(self.model.variables['12x!!@#5_3'].ub, 666) + repickled = pickle.loads(pickle.dumps(self.model)) + var_from_pickle = repickled.variables['12x!!@#5_3'] + self.assertEqual(var_from_pickle.name, repickled.problem.getVarByName(var.name).VarName) + + def test_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') + self.assertEqual(var.problem, self.model) + self.model.remove(var) + self.model.problem.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') + self.assertNotIn(var, self.model.variables.values()) + self.assertEqual(self.model.problem.getVarByName(var.name), None) + self.assertEqual(var.problem, None) + + def test_remove_variable_str(self): + var = self.model.variables.values()[0] + self.model.remove(var.name) + self.model.problem.update() + self.assertNotIn(var, self.model.variables.values()) + self.assertEqual(self.model.problem.getVarByName(var.name), None) + self.assertEqual(var.problem, None) + + def test_add_constraints(self): + x = Variable('x', lb=0, ub=1, type='binary') + y = Variable('y', lb=-181133.3, ub=12000., type='continuous') + z = Variable('z', lb=0., ub=10., type='integer') + constr1 = Constraint(0.3 * x + 0.4 * y + 66. * z, lb=-100, ub=0., name='test') + constr2 = Constraint(2.333 * x + y + 3.333, ub=100.33, name='test2') + constr3 = Constraint(2.333 * x + y + z, lb=-300) + constr4 = Constraint(x, lb=-300, ub=-300) + constr5 = Constraint(3*x) + self.model.add(constr1) + self.model.add(constr2) + self.model.add(constr3) + self.model.add(constr4) + self.model.problem.update() + self.assertIn(constr1.name, self.model.constraints) + self.assertIn(constr2.name, self.model.constraints) + self.assertIn(constr3.name, self.model.constraints) + self.assertIn(constr4.name, self.model.constraints) + # constr1 + coeff_dict = dict() + internal_constraint = self.model.problem.getConstrByName(constr1.name) + row = self.model.problem.getRow(internal_constraint) + for i in range(row.size()): + coeff_dict[row.getVar(i).VarName] = row.getCoeff(i) + self.assertDictEqual(coeff_dict, {'x': 0.3, 'y': 0.4, 'z': 66., 'test_aux': -1.0}) + self.assertEqual(internal_constraint.RHS, 0) + self.assertEqual(self.model.problem.getVarByName(internal_constraint.getAttr('ConstrName')+'_aux'), 100) + # constr2 + coeff_dict = dict() + internal_constraint = self.model.problem.getConstrByName(constr2.name) + row = self.model.problem.getRow(internal_constraint) + for i in range(row.size()): + coeff_dict[row.getVar(i).VarName] = row.getCoeff(i) + self.assertDictEqual(coeff_dict, {'x': 2.333, 'y': 1.}) + self.assertEqual(internal_constraint.RHS, 96.997) + self.assertEqual(internal_constraint.Sense, '<') + # constr3 + coeff_dict = dict() + internal_constraint = self.model.problem.getConstrByName(constr3.name) + print(internal_constraint) + row = self.model.problem.getRow(internal_constraint) + for i in range(row.size()): + coeff_dict[row.getVar(i).VarName] = row.getCoeff(i) + self.assertDictEqual(coeff_dict, {'x': 2.333, 'y': 1., 'z': 1.}) + self.assertEqual(internal_constraint.RHS, -300) + self.assertEqual(internal_constraint.Sense, '>') + # constr4 + coeff_dict = dict() + internal_constraint = self.model.problem.getConstrByName(constr4.name) + print(internal_constraint) + row = self.model.problem.getRow(internal_constraint) + for i in range(row.size()): + coeff_dict[row.getVar(i).VarName] = row.getCoeff(i) + self.assertDictEqual(coeff_dict, {'x': 1}) + self.assertEqual(internal_constraint.RHS, -300) + self.assertEqual(internal_constraint.Sense, '=') + + def test_remove_constraints(self): + x = Variable('x', type='binary') + y = Variable('y', lb=-181133.3, ub=12000., type='continuous') + z = Variable('z', lb=3, ub=3, type='integer') + constr1 = Constraint(0.3 * x + 0.4 * y + 66. * z, lb=-100, ub=0., name='test') + self.assertEqual(constr1.problem, None) + self.model.add(constr1) + self.model.update() + self.assertEqual(constr1.problem, self.model) + self.assertIn(constr1.name, self.model.constraints) + print('test', constr1.name in self.model.constraints.keys()) + self.model.remove(constr1.name) + self.model.update() + self.assertEqual(constr1.problem, None) + self.assertNotIn(constr1, self.model.constraints) + + def test_add_nonlinear_constraint_raises(self): + x = Variable('x', type='binary') + y = Variable('y', lb=-181133.3, ub=12000., type='continuous') + z = Variable('z', lb=10, type='integer') + c = Constraint(0.3 * x + 0.4 * y ** 2 + 66. * z, lb = -100, ub = 0., name = 'test') + self.model.add(c) + self.assertRaisesRegexp(ValueError, 'GUROBI currently only supports linear constraint.*', self.model.update) + + def test_change_of_constraint_is_reflected_in_low_level_solver(self): + x = Variable('x', lb=-83.3, ub=1324422.) + y = Variable('y', lb=-181133.3, ub=12000.) + 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') + 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') + + 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.model.update() + grb_obj = self.model.problem.getObjective() + grb_x = self.model.problem.getVarByName(x.name) + grb_y = self.model.problem.getVarByName(y.name) + expected = {grb_x: 0.3, grb_y: 0.4} + for i in range(grb_obj.size()): + self.assertEqual(grb_obj.getCoeff(i), expected[grb_obj.getVar(i)]) + z = Variable('z', lb=4, ub=4, type='integer') + grb_z = self.model.problem.getVarByName(z.name) + self.model.objective += 77. * z + expected[grb_z] = 77. + self.model.update() + for i in range(grb_obj.size()): + self.assertEqual(grb_obj.getCoeff(i), expected[grb_obj.getVar(i)]) + + def test_change_variable_bounds(self): + inner_prob = self.model.problem + inner_problem_bounds = [(variable.LB, variable.UB) for variable in inner_prob.getVars()] + bounds = [(var.lb, var.ub) for var in self.model.variables.values()] + self.assertEqual(bounds, inner_problem_bounds) + for var in self.model.variables.values(): + var.lb = random.uniform(-1000, 1000) + var.ub = random.uniform(var.lb, 1000) + self.model.update() + inner_problem_bounds_new = [(variable.LB, variable.UB) for variable in inner_prob.getVars()] + bounds_new = [(var.lb, var.ub) for var in self.model.variables.values()] + self.assertNotEqual(bounds, bounds_new) + self.assertNotEqual(inner_problem_bounds, inner_problem_bounds_new) + self.assertEqual(bounds_new, inner_problem_bounds_new) + + def test_change_variable_type(self): + for variable in self.model.variables: + variable.type = 'integer' + self.model.update() + for variable in self.model.problem.getVars(): + self.assertEqual(variable.VType, gurobipy.GRB.INTEGER) + + def test_change_constraint_bounds(self): + inner_prob = self.model.problem + inner_problem_bounds = [] + for constr in inner_prob.getConstrs(): + aux_var = inner_prob.getVarByName(constr.getAttr('ConstrName') + '_aux') + if aux_var is None: + inner_problem_bounds.append((constr.RHS, constr.RHS)) + else: + inner_problem_bounds.append((aux_var.UB, constr.RHS)) + print(len(self.model.constraints)) + print(len(self.model.problem.getConstrs())) + bounds = [(constr.lb, constr.ub) for constr in self.model.constraints] + print('bounds', inner_problem_bounds) + 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') + + def test_optimize(self): + self.model.optimize() + self.assertEqual(self.model.status, 'optimal') + self.assertAlmostEqual(self.model.objective.value, 0.8739215069684303) + + def test_optimize_milp(self): + problem = gurobipy.read(TESTMILPMODELPATH) + milp_model = Model(problem=problem) + milp_model.optimize() + self.assertEqual(milp_model.status, 'optimal') + self.assertAlmostEqual(milp_model.objective.value, 122.5) + for variable in milp_model.variables: + if variable.type == 'integer': + self.assertEqual(variable.primal % 1, 0) + + def test_change_objective(self): + """Test that all different kinds of linear objective specification work.""" + print(self.model.variables.values()[0:2]) + v1, v2 = self.model.variables.values()[0:2] + self.model.objective = Objective(1. * v1 + 1. * v2) + self.assertEqual(self.model.objective.__str__(), 'Maximize\n1.0*R_PGK + 1.0*R_Biomass_Ecoli_core_w_GAM') + self.model.objective = Objective(v1 + v2) + self.assertEqual(self.model.objective.__str__(), 'Maximize\n1.0*R_PGK + 1.0*R_Biomass_Ecoli_core_w_GAM') + + @unittest.skip('Incomplete') + def test_number_objective(self): + self.model.objective = Objective(0.) + self.model.update() + self.assertEqual(self.model.objective.__str__(), 'Maximize\n0') + obj_coeff = list() + print(self.model.problem.getObjective().size()) + print(self.model.problem.getObjective().getVar(0)) + # for i in range(1, glp_get_num_cols(self.model.problem) + 1): + # obj_coeff.append(glp_get_obj_coef(self.model.problem, i)) + self.assertEqual(set(obj_coeff), {0.}) + + @unittest.skip('Incomplete') + def test_raise_on_non_linear_objective(self): + """Test that an exception is raised when a non-linear objective is added to the model.""" + v1, v2 = self.model.variables.values()[0:2] + self.assertRaises(ValueError, Objective, v1*v2) + + @unittest.skip('Not supported yet') + def test_iadd_objective(self): + v2, v3 = self.model.variables.values()[1:3] + print(v2, v3) + # 1/0 + self.model.objective += 2. * v2 - 3. * v3 + internal_objective = self.model.problem.getObjective() + result = {} + for i in range(internal_objective.size()): + var = internal_objective.getVar(i) + coeff = internal_objective.getCoeff(i) + result[var.VarName] = coeff + self.assertDictEqual(result, {'R_Biomass_Ecoli_core_w_GAM': 1.0}) + self.model.update() + self.assertDictEqual(result, {'R_Biomass_Ecoli_core_w_GAM': 1.0, 'R_PGK': 2, 'R_GAPD': -3}) + + @unittest.skip('Not supported yet') + def test_imul_objective(self): + self.model.objective *= 2. + obj_coeff = list() + for i in range(len(self.model.variables)): + obj_coeff.append(glp_get_obj_coef(self.model.problem, i)) + self.assertEqual(obj_coeff, + [0.0, 2.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.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.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.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.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0] + ) + + 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') + + def test_timeout(self): + self.model.configuration.timeout = 0 + status = self.model.optimize() + print(status) + self.assertEqual(status, 'time_limit') + + def test_instantiating_model_with_non_glpk_problem_raises(self): + self.assertRaises(TypeError, Model, problem='Chicken soup') + + def test_primal_values(self): + self.model.optimize() + for k, v in self.model.primal_values.items(): + self.assertEquals(v, self.model.variables[k].primal) + + def test_reduced_costs(self): + self.model.optimize() + for k, v in self.model.reduced_costs.items(): + self.assertEquals(v, self.model.variables[k].dual) + + def test_dual_values(self): + self.model.optimize() + for k, v in self.model.dual_values.items(): + self.assertEquals(v, self.model.constraints[k].primal) + + def test_shadow_prices(self): + self.model.optimize() + for k, v in self.model.shadow_prices.items(): + self.assertEquals(v, self.model.constraints[k].dual) + + def test_clone_solver(self): + self.assertEquals(self.model.configuration.verbosity, 0) + self.model.configuration.verbosity = 3 + cloned_model = Model.clone(self.model) + self.assertEquals(cloned_model.configuration.verbosity, 3) + self.assertEquals(len(cloned_model.variables), len(self.model.variables)) + self.assertEquals(len(cloned_model.constraints), len(self.model.constraints)) + + +if __name__ == '__main__': + nose.runmodule() diff --git a/tests/test_netlib_gurobi_interface.py b/tests/test_netlib_gurobi_interface.py new file mode 100644 index 00000000..c5d1e036 --- /dev/null +++ b/tests/test_netlib_gurobi_interface.py @@ -0,0 +1,168 @@ +# Copyright (c) 2013 Novo Nordisk Foundation Center for Biosustainability, DTU. +# See LICENSE for details. + +import unittest +import pickle +import tempfile +import glob +import tarfile +from functools import partial + +import os +import nose + +import six + +import logging +log = logging.getLogger(__name__) +log.setLevel(logging.DEBUG) + +if six.PY3: + nose.SkipTest('Skipping because py3') +else: + try: + import gurobipy + + from optlang.gurobi_interface import Model + + + with open(os.path.join(os.path.dirname(__file__), 'data/the_final_netlib_results.pcl'), 'rb') as fhandle: + THE_FINAL_NETLIB_RESULTS = pickle.load(fhandle) + + + class TestClass(object): + """docstring for TestClass""" + + def __init__(self, arg): + super(TestClass, self).__init__() + self.arg = arg + + + # noinspection PyShadowingNames + def read_netlib_sif_gurobi(fhandle): + tmp_file = tempfile.mktemp(suffix='.mps') + with open(tmp_file, 'w') as tmp_handle: + content = ''.join([str(s) for s in fhandle if str(s).strip()]) + tmp_handle.write(content) + fhandle.close() + log.debug(tmp_file) + problem = gurobipy.read(tmp_file) + return problem + + + def check_dimensions(problem, model): + """ + Tests that the glpk problem and the interface model have the same + number of rows (constraints) and columns (variables). + """ + assert len(problem.getVars()) == len(model.variables) + + + def check_objval(problem, model_objval): + """ + Check that ... + """ + if problem.getAttr('Status') == gurobipy.GRB.OPTIMAL: + problem_objval = problem.getAttr('ObjVal') + else: + problem_objval = None + nose.tools.assert_almost_equal(problem_objval, model_objval, places=4) + + + def check_objval_against_the_final_netlib_results(netlib_id, model_objval): + relative_error = abs(1 - (model_objval / float(THE_FINAL_NETLIB_RESULTS[netlib_id]['Objvalue']))) + print(relative_error) + nose.tools.assert_true(relative_error < 0.01) + # nose.tools.assert_almost_equal(model_objval, float(THE_FINAL_NETLIB_RESULTS[netlib_id]['Objvalue']), places=4) + + + def test_netlib(netlib_tar_path=os.path.join(os.path.dirname(__file__), 'data/netlib_lp_problems.tar.gz')): + """ + Test netlib with glpk interface + """ + tar = tarfile.open(netlib_tar_path) + model_paths_in_tar = glob.fnmatch.filter(tar.getnames(), '*.SIF') + + for model_path_in_tar in model_paths_in_tar: + print(model_path_in_tar) + netlib_id = os.path.basename(model_path_in_tar).replace('.SIF', '') + # TODO: get the following problems to work + # E226 seems to be a MPS related problem, see http://lists.gnu.org/archive/html/bug-glpk/2003-01/msg00003.html + if netlib_id in ('AGG', 'E226', 'SCSD6', 'BLEND', 'DFL001', 'FORPLAN', 'GFRD-PNC', 'SIERRA'): + # def test_skip(netlib_id): + # raise SkipTest('Skipping netlib problem %s ...' % netlib_id) + # test_skip(netlib_id) + # class TestWeirdNetlibProblems(unittest.TestCase): + + # @unittest.skip('Skipping netlib problem') + # def test_fail(): + # pass + continue + # TODO: For now, test only models that are covered by the final netlib results + else: + if netlib_id not in THE_FINAL_NETLIB_RESULTS.keys(): + continue + fhandle = tar.extractfile(model_path_in_tar) + problem = read_netlib_sif_gurobi(fhandle) + model = Model(problem=problem) + model.configuration.presolve = True + model.configuration.verbosity = 3 + func = partial(check_dimensions, problem, model) + func.description = "test_netlib_check_dimensions_%s (%s)" % (netlib_id, os.path.basename(str(__file__))) + yield func + + model.optimize() + if model.status == 'optimal': + model_objval = model.objective.value + else: + raise Exception('No optimal solution found for netlib model %s' % netlib_id) + + func = partial(check_objval, problem, model_objval) + func.description = "test_netlib_check_objective_value_%s (%s)" % ( + netlib_id, os.path.basename(str(__file__))) + yield func + + func = partial(check_objval_against_the_final_netlib_results, netlib_id, model_objval) + func.description = "test_netlib_check_objective_value__against_the_final_netlib_results_%s (%s)" % ( + netlib_id, os.path.basename(str(__file__))) + yield func + + except ImportError as e: + + if str(e).find('gurobipy') >= 0: + class TestMissingDependency(unittest.TestCase): + + @unittest.skip('Missing dependency - ' + str(e)) + def test_fail(self): + pass + else: + raise + +if __name__ == '__main__': + # tar = tarfile.open('data/netlib_lp_problems.tar.gz') + # model_paths_in_tar = glob.fnmatch.filter(tar.getnames(), '*.SIF') + # for model_path_in_tar in model_paths_in_tar: + # try: + # + # fhandle = tar.extractfile(model_path_in_tar) + # problem = read_netlib_sif_cplex(fhandle) + # + # except cplex.exceptions.CplexSolverError, e: + # print model_path_in_tar + # print problem + # print problem.get_problem_name() + # print e + # fhandle = tar.extractfile('./netlib/ADLITTLE.SIF') + # glpk_problem = read_netlib_sif_glpk(fhandle) + # glp_simplex(glpk_problem, None) + # print glp_get_obj_val(glpk_problem) + # print glpk_problem + # fhandle = tar.extractfile('./netlib/ADLITTLE.SIF') + # glpk_problem = read_netlib_sif_glpk(fhandle) + # model = Model(problem=glpk_problem) + # glp_simplex(glpk_problem, None) + # model.optimize() + # print model.objective.value + # print model + # test_netlib().next() + nose.runmodule() \ No newline at end of file