Merge pull request #791 from opencobra/fix-inf-bounds

Fix inf bounds

cobra/core/model.py

@@ -16,7 +16,7 @@
 from cobra.core.configuration import Configuration
 from cobra.core.dictlist import DictList
 from cobra.core.object import Object
-from cobra.core.reaction import Reaction, separate_forward_and_reverse_bounds
+from cobra.core.reaction import Reaction
 from cobra.core.solution import get_solution
 from cobra.exceptions import SolverNotFound
 from cobra.medium import find_boundary_types, sbo_terms
@@ -806,25 +806,14 @@ def _populate_solver(self, reaction_list, metabolite_list=None):
         self.add_cons_vars(to_add)
 
         for reaction in reaction_list:
-
             if reaction.id not in self.variables:
-
-                reverse_lb, reverse_ub, forward_lb, forward_ub = \
-                    separate_forward_and_reverse_bounds(*reaction.bounds)
-
-                forward_variable = self.problem.Variable(
-                    reaction.id, lb=forward_lb, ub=forward_ub)
-                reverse_variable = self.problem.Variable(
-                    reaction.reverse_id, lb=reverse_lb, ub=reverse_ub)
-
+                forward_variable = self.problem.Variable(reaction.id)
+                reverse_variable = self.problem.Variable(reaction.reverse_id)
                 self.add_cons_vars([forward_variable, reverse_variable])
-
             else:
-
                 reaction = self.reactions.get_by_id(reaction.id)
                 forward_variable = reaction.forward_variable
                 reverse_variable = reaction.reverse_variable
-
             for metabolite, coeff in six.iteritems(reaction.metabolites):
                 if metabolite.id in self.constraints:
                     constraint = self.constraints[metabolite.id]
@@ -834,11 +823,13 @@ def _populate_solver(self, reaction_list, metabolite_list=None):
                         name=metabolite.id,
                         lb=0, ub=0)
                     self.add_cons_vars(constraint, sloppy=True)
-
                 constraint_terms[constraint][forward_variable] = coeff
                 constraint_terms[constraint][reverse_variable] = -coeff
 
         self.solver.update()
+        for reaction in reaction_list:
+            reaction = self.reactions.get_by_id(reaction.id)
+            reaction.update_variable_bounds()
         for constraint, terms in six.iteritems(constraint_terms):
             constraint.set_linear_coefficients(terms)
 

cobra/core/reaction.py

@@ -7,6 +7,7 @@
 from collections import defaultdict
 from copy import copy, deepcopy
 from functools import partial
+from math import isinf
 from operator import attrgetter
 from warnings import warn
 
@@ -91,9 +92,6 @@ def __init__(self, id=None, name='', subsystem='', lower_bound=0.0,
         self._upper_bound = upper_bound if upper_bound is not None else \
             CONFIGURATION.upper_bound
 
-        self._reverse_variable = None
-        self._forward_variable = None
-
     def _set_id_with_model(self, value):
         if value in self.model.reactions:
             raise ValueError("The model already contains a reaction with"
@@ -186,6 +184,39 @@ def __deepcopy__(self, memo):
         cop = deepcopy(super(Reaction, self), memo)
         return cop
 
+    @staticmethod
+    def _check_bounds(lb, ub):
+        if lb > ub:
+            raise ValueError(
+                "The lower bound must be less than or equal to the upper "
+                "bound ({} <= {}).".format(lb, ub))
+
+    def update_variable_bounds(self):
+        if self.model is None:
+            return
+        # We know that `lb <= ub`.
+        if self._lower_bound > 0:
+            self.forward_variable.set_bounds(
+                lb=None if isinf(self._lower_bound) else self._lower_bound,
+                ub=None if isinf(self._upper_bound) else self._upper_bound
+            )
+            self.reverse_variable.set_bounds(lb=0, ub=0)
+        elif self._upper_bound < 0:
+            self.forward_variable.set_bounds(lb=0, ub=0)
+            self.reverse_variable.set_bounds(
+                lb=None if isinf(self._upper_bound) else -self._upper_bound,
+                ub=None if isinf(self._lower_bound) else -self._lower_bound
+            )
+        else:
+            self.forward_variable.set_bounds(
+                lb=0,
+                ub=None if isinf(self._upper_bound) else self._upper_bound
+            )
+            self.reverse_variable.set_bounds(
+                lb=0,
+                ub=None if isinf(self._lower_bound) else -self._lower_bound
+            )
+
     @property
     def lower_bound(self):
         """Get or set the lower bound
@@ -203,11 +234,18 @@ def lower_bound(self):
     @lower_bound.setter
     @resettable
     def lower_bound(self, value):
-        if self.upper_bound < value:
+        if self._upper_bound < value:
+            warn("You are constraining the reaction '{}' to a fixed flux "
+                 "value of {}. Did you intend to do this? We are planning to "
+                 "remove this behavior in a future release. Please let us "
+                 "know your opinion at "
+                 "https://github.com/opencobra/cobrapy/issues/793."
+                 "".format(self.id, value), DeprecationWarning)
             self.upper_bound = value
-
+        # Validate bounds before setting them.
+        self._check_bounds(value, self._upper_bound)
         self._lower_bound = value
-        update_forward_and_reverse_bounds(self, 'lower')
+        self.update_variable_bounds()
 
     @property
     def upper_bound(self):
@@ -226,11 +264,18 @@ def upper_bound(self):
     @upper_bound.setter
     @resettable
     def upper_bound(self, value):
-        if self.lower_bound > value:
+        if self._lower_bound > value:
+            warn("You are constraining the reaction '{}' to a fixed flux "
+                 "value of {}. Did you intend to do this? We are planning to "
+                 "remove this behavior in a future release. Please let us "
+                 "know your opinion at "
+                 "https://github.com/opencobra/cobrapy/issues/793."
+                 "".format(self.id, value), DeprecationWarning)
             self.lower_bound = value
-
+        # Validate bounds before setting them.
+        self._check_bounds(self._lower_bound, value)
         self._upper_bound = value
-        update_forward_and_reverse_bounds(self, 'upper')
+        self.update_variable_bounds()
 
     @property
     def bounds(self):
@@ -248,9 +293,12 @@ def bounds(self):
     @bounds.setter
     @resettable
     def bounds(self, value):
-        assert value[0] <= value[1], "Invalid bounds: {}".format(value)
-        self._lower_bound, self._upper_bound = value
-        update_forward_and_reverse_bounds(self)
+        lower, upper = value
+        # Validate bounds before setting them.
+        self._check_bounds(lower, upper)
+        self._lower_bound = lower
+        self._upper_bound = upper
+        self.update_variable_bounds()
 
     @property
     def flux(self):
@@ -1083,64 +1131,3 @@ def _repr_html_(self):
                        self.build_reaction_string(True), 200),
                    gpr=format_long_string(self.gene_reaction_rule, 100),
                    lb=self.lower_bound, ub=self.upper_bound)
-
-
-def separate_forward_and_reverse_bounds(lower_bound, upper_bound):
-    """Split a given (lower_bound, upper_bound) interval into a negative
-    component and a positive component. Negative components are negated
-    (returns positive ranges) and flipped for usage with forward and reverse
-    reactions bounds
-
-    Parameters
-    ----------
-    lower_bound : float
-        The lower flux bound
-    upper_bound : float
-        The upper flux bound
-    """
-
-    assert lower_bound <= upper_bound, "lower bound is greater than upper"
-
-    bounds_list = [0, 0, lower_bound, upper_bound]
-    bounds_list.sort()
-
-    return -bounds_list[1], -bounds_list[0], bounds_list[2], bounds_list[3]
-
-
-def update_forward_and_reverse_bounds(reaction, direction='both'):
-    """For the given reaction, update the bounds in the forward and
-    reverse variable bounds.
-
-    Parameters
-    ----------
-    reaction : cobra.Reaction
-       The reaction to operate on
-    direction : string
-       Either 'both', 'upper' or 'lower' for updating the corresponding flux
-       bounds.
-    """
-
-    reverse_lb, reverse_ub, forward_lb, forward_ub = \
-        separate_forward_and_reverse_bounds(*reaction.bounds)
-
-    try:
-        # Clear the original bounds to avoid complaints
-        if direction == 'both':
-            reaction.forward_variable._ub = None
-            reaction.reverse_variable._lb = None
-            reaction.reverse_variable._ub = None
-            reaction.forward_variable._lb = None
-
-            reaction.forward_variable.set_bounds(lb=forward_lb, ub=forward_ub)
-            reaction.reverse_variable.set_bounds(lb=reverse_lb, ub=reverse_ub)
-
-        elif direction == 'upper':
-            reaction.forward_variable.ub = forward_ub
-            reaction.reverse_variable.lb = reverse_lb
-
-        elif direction == 'lower':
-            reaction.reverse_variable.ub = reverse_ub
-            reaction.forward_variable.lb = forward_lb
-
-    except AttributeError:
-        pass

cobra/manipulation/__init__.py

@@ -7,8 +7,6 @@
     delete_model_genes, find_gene_knockout_reactions, remove_genes,
     undelete_model_genes)
 from cobra.manipulation.modify import (
-    convert_to_irreversible, escape_ID, get_compiled_gene_reaction_rules,
-    revert_to_reversible)
+    escape_ID, get_compiled_gene_reaction_rules)
 from cobra.manipulation.validate import (
-    check_mass_balance, check_metabolite_compartment_formula,
-    check_reaction_bounds)
+    check_mass_balance, check_metabolite_compartment_formula)

cobra/manipulation/delete.py

@@ -6,7 +6,6 @@
 
 from six import iteritems, string_types
 
-from cobra.core import Model
 from cobra.core.gene import ast2str, eval_gpr, parse_gpr
 
 

cobra/manipulation/modify.py

@@ -8,7 +8,7 @@
 
 from six import iteritems
 
-from cobra.core import Gene, Reaction
+from cobra.core import Reaction
 from cobra.core.gene import ast2str
 from cobra.manipulation.delete import get_compiled_gene_reaction_rules
 from cobra.util.solver import set_objective
@@ -116,78 +116,3 @@ def visit_Name(self, node):
         rxn.gene_reaction_rule = rxn._gene_reaction_rule
     for i in remove_genes:
         cobra_model.genes.remove(i)
-
-
-def convert_to_irreversible(cobra_model):
-    """Split reversible reactions into two irreversible reactions
-
-    These two reactions will proceed in opposite directions. This
-    guarentees that all reactions in the model will only allow
-    positive flux values, which is useful for some modeling problems.
-
-    cobra_model: A Model object which will be modified in place.
-
-    """
-    warn("deprecated, not applicable for optlang solvers", DeprecationWarning)
-    reactions_to_add = []
-    coefficients = {}
-    for reaction in cobra_model.reactions:
-        # If a reaction is reverse only, the forward reaction (which
-        # will be constrained to 0) will be left in the model.
-        if reaction.lower_bound < 0:
-            reverse_reaction = Reaction(reaction.id + "_reverse")
-            reverse_reaction.lower_bound = max(0, -reaction.upper_bound)
-            reverse_reaction.upper_bound = -reaction.lower_bound
-            coefficients[
-                reverse_reaction] = reaction.objective_coefficient * -1
-            reaction.lower_bound = max(0, reaction.lower_bound)
-            reaction.upper_bound = max(0, reaction.upper_bound)
-            # Make the directions aware of each other
-            reaction.notes["reflection"] = reverse_reaction.id
-            reverse_reaction.notes["reflection"] = reaction.id
-            reaction_dict = {k: v * -1
-                             for k, v in iteritems(reaction._metabolites)}
-            reverse_reaction.add_metabolites(reaction_dict)
-            reverse_reaction._model = reaction._model
-            reverse_reaction._genes = reaction._genes
-            for gene in reaction._genes:
-                gene._reaction.add(reverse_reaction)
-            reverse_reaction.subsystem = reaction.subsystem
-            reverse_reaction._gene_reaction_rule = reaction._gene_reaction_rule
-            reactions_to_add.append(reverse_reaction)
-    cobra_model.add_reactions(reactions_to_add)
-    set_objective(cobra_model, coefficients, additive=True)
-
-
-def revert_to_reversible(cobra_model, update_solution=True):
-    """This function will convert an irreversible model made by
-    convert_to_irreversible into a reversible model.
-
-    cobra_model : cobra.Model
-        A model which will be modified in place.
-    update_solution: bool
-        This option is ignored since `model.solution` was removed.
-    """
-    warn("deprecated, not applicable for optlang solvers", DeprecationWarning)
-    reverse_reactions = [x for x in cobra_model.reactions
-                         if "reflection" in x.notes and
-                         x.id.endswith('_reverse')]
-
-    # If there are no reverse reactions, then there is nothing to do
-    if len(reverse_reactions) == 0:
-        return
-
-    for reverse in reverse_reactions:
-        forward_id = reverse.notes.pop("reflection")
-        forward = cobra_model.reactions.get_by_id(forward_id)
-        forward.lower_bound = -reverse.upper_bound
-        if forward.upper_bound == 0:
-            forward.upper_bound = -reverse.lower_bound
-
-        if "reflection" in forward.notes:
-            forward.notes.pop("reflection")
-
-    # Since the metabolites and genes are all still in
-    # use we can do this faster removal step.  We can
-    # probably speed things up here.
-    cobra_model.remove_reactions(reverse_reactions)

cobra/manipulation/validate.py

@@ -2,9 +2,6 @@
 
 from __future__ import absolute_import
 
-from math import isinf, isnan
-from warnings import warn
-
 
 NOT_MASS_BALANCED_TERMS = {"SBO:0000627",  # EXCHANGE
                            "SBO:0000628",  # DEMAND
@@ -24,30 +21,6 @@ def check_mass_balance(model):
     return unbalanced
 
 
-# no longer strictly necessary, done by optlang solver interfaces
-def check_reaction_bounds(model):
-    warn("no longer necessary, done by optlang solver interfaces",
-         DeprecationWarning)
-    errors = []
-    for reaction in model.reactions:
-        if reaction.lower_bound > reaction.upper_bound:
-            errors.append("Reaction '%s' has lower bound > upper bound" %
-                          reaction.id)
-        if isinf(reaction.lower_bound):
-            errors.append("Reaction '%s' has infinite lower_bound" %
-                          reaction.id)
-        elif isnan(reaction.lower_bound):
-            errors.append("Reaction '%s' has NaN for lower_bound" %
-                          reaction.id)
-        if isinf(reaction.upper_bound):
-            errors.append("Reaction '%s' has infinite upper_bound" %
-                          reaction.id)
-        elif isnan(reaction.upper_bound):
-            errors.append("Reaction '%s' has NaN for upper_bound" %
-                          reaction.id)
-    return errors
-
-
 def check_metabolite_compartment_formula(model):
     errors = []
     for met in model.metabolites:

cobra/test/test_core/test_core_reaction.py

@@ -272,7 +272,7 @@ def test_build_from_string(model):
 
 def test_bounds_setter(model):
     rxn = model.reactions.get_by_id("PGI")
-    with pytest.raises(AssertionError):
+    with pytest.raises(ValueError):
         rxn.bounds = (1, 0)