cobra/__init__.py

@@ -13,7 +13,7 @@ from cobra.core import (
     DictList, Gene, Metabolite, Model, Object, Reaction, Species)
 from cobra.util.version_info import show_versions
 
-__version__ = "0.11.0"
+__version__ = "0.11.3"
 
 # set the warning format to be prettier and fit on one line
 _cobra_path = _dirname(_abspath(__file__))

cobra/core/reaction.py

@@ -420,7 +420,7 @@ class Reaction(Object):
             if g not in self._genes:  # if an old gene is not a new gene
                 try:
                     g._reaction.remove(self)
-                except:
+                except KeyError:
                     warn("could not remove old gene %s from reaction %s" %
                          (g.id, self.id))
 
@@ -618,6 +618,7 @@ class Reaction(Object):
     __radd__ = __add__
 
     def __iadd__(self, other):
+
         self.add_metabolites(other._metabolites, combine=True)
         gpr1 = self.gene_reaction_rule.strip()
         gpr2 = other.gene_reaction_rule.strip()
@@ -643,11 +644,26 @@ class Reaction(Object):
     def __imul__(self, coefficient):
         """Scale coefficients in a reaction by a given value
 
-        E.g. A -> B becomes 2A -> 2B
+        E.g. A -> B becomes 2A -> 2B.
+
+        If coefficient is less than zero, the reaction is reversed and the
+        bounds are swapped.
         """
         self._metabolites = {
             met: value * coefficient
             for met, value in iteritems(self._metabolites)}
+
+        if coefficient < 0:
+            self.bounds = (-self.upper_bound, -self.lower_bound)
+
+        if self._model:
+            self._model._populate_solver([self])
+
+        context = get_context(self)
+        if context:
+            context(partial(self._model._populate_solver, [self]))
+            context(partial(self.__imul__, 1./coefficient))
+
         return self
 
     def __mul__(self, coefficient):
@@ -757,37 +773,25 @@ class Reaction(Object):
                 # reaction
                 metabolite._reaction.add(self)
 
-        for metabolite, the_coefficient in list(self._metabolites.items()):
-            if the_coefficient == 0:
-                # make the metabolite aware that it no longer participates
-                # in this reaction
-                metabolite._reaction.remove(self)
-                self._metabolites.pop(metabolite)
-
         # from cameo ...
         model = self.model
         if model is not None:
             model.add_metabolites(new_metabolites)
 
-            for metabolite, coefficient in metabolites_to_add.items():
-
-                if isinstance(metabolite,
-                              str):  # support metabolites added as strings.
-                    metabolite = model.metabolites.get_by_id(metabolite)
-                if combine:
-                    try:
-                        old_coefficient = old_coefficients[metabolite]
-                    except KeyError:
-                        pass
-                    else:
-                        coefficient = coefficient + old_coefficient
-
+            for metabolite, coefficient in self._metabolites.items():
                 model.constraints[
                     metabolite.id].set_linear_coefficients(
                     {self.forward_variable: coefficient,
                      self.reverse_variable: -coefficient
                      })
 
+        for metabolite, the_coefficient in list(self._metabolites.items()):
+            if the_coefficient == 0:
+                # make the metabolite aware that it no longer participates
+                # in this reaction
+                metabolite._reaction.remove(self)
+                self._metabolites.pop(metabolite)
+
         context = get_context(self)
         if context and reversibly:
             if combine:

cobra/flux_analysis/parsimonious.py

@@ -94,7 +94,7 @@ def add_pfba(model, objective=None, fraction_of_optimum=1.0):
     if objective is not None:
         model.objective = objective
     if model.solver.objective.name == '_pfba_objective':
-        raise ValueError('model already has pfba objective')
+        raise ValueError('The model already has a pFBA objective.')
     sutil.fix_objective_as_constraint(model, fraction=fraction_of_optimum)
     reaction_variables = ((rxn.forward_variable, rxn.reverse_variable)
                           for rxn in model.reactions)

cobra/flux_analysis/sampling.py

@@ -24,18 +24,14 @@ from cobra.util import (create_stoichiometric_matrix, constraint_matrices,
 LOGGER = getLogger(__name__)
 """The logger for the package."""
 
-bounds_tol = np.finfo(np.float32).eps
+bounds_tol = 1e-6
 """The tolerance used for checking bounds feasibility."""
 
-feasibility_tol = bounds_tol
+feasibility_tol = 1e-6
 """The tolerance used for checking equalities feasibility."""
 
-nproj = 1000000
-"""Reproject the solution into the feasibility space every nproj iterations."""
-
-nproj_center = 10000
-"""Reproject the center into the nullspace every nproj_center iterations.
-   Only used for inhomogeneous problems."""
+max_tries = 100
+"""Maximum number of retries for sampling."""
 
 Problem = namedtuple("Problem",
                      ["equalities", "b", "inequalities", "bounds",
@@ -106,7 +102,7 @@ def shared_np_array(shape, data=None, integer=False):
 
 
 # Has to be declared outside of class to be used for multiprocessing :(
-def _step(sampler, x, delta, fraction=None):
+def _step(sampler, x, delta, fraction=None, tries=0):
     """Sample a new feasible point from the point `x` in direction `delta`."""
     prob = sampler.problem
     valid = ((np.abs(delta) > feasibility_tol) &
@@ -116,14 +112,19 @@ def _step(sampler, x, delta, fraction=None):
     valphas = (valphas / delta[valid]).flatten()
     if prob.bounds.shape[0] > 0:
         # permissible alphas for staying in constraint bounds
+        ineqs = prob.inequalities.dot(delta)
+        valid = np.abs(ineqs) > feasibility_tol
         balphas = ((1.0 - bounds_tol) * prob.bounds -
-                   prob.inequalities.dot(x))
-        balphas = (balphas / prob.inequalities.dot(delta)).flatten()
+                   prob.inequalities.dot(x))[:, valid]
+        balphas = (balphas / ineqs[valid]).flatten()
         # combined alphas
         alphas = np.hstack([valphas, balphas])
     else:
         alphas = valphas
-    alpha_range = (alphas[alphas > 0.0].min(), alphas[alphas <= 0.0].max())
+    pos_alphas = alphas[alphas > 0.0]
+    neg_alphas = alphas[alphas <= 0.0]
+    alpha_range = np.array([neg_alphas.max() if len(neg_alphas) > 0 else 0,
+                            pos_alphas.min() if len(pos_alphas) > 0 else 0])
 
     if fraction:
         alpha = alpha_range[0] + fraction * (alpha_range[1] - alpha_range[0])
@@ -133,12 +134,21 @@ def _step(sampler, x, delta, fraction=None):
 
     # Numerical instabilities may cause bounds invalidation
     # reset sampler and sample from one of the original warmup directions
-    # if that occurs
-    if np.any(sampler._bounds_dist(p) < -bounds_tol):
+    # if that occurs. Also reset if we got stuck.
+    if (np.any(sampler._bounds_dist(p) < -bounds_tol) or
+            np.abs(np.abs(alpha_range).max() * delta).max() < bounds_tol):
+        if tries > max_tries:
+            raise RuntimeError("Can not escape sampling region, model seems"
+                               " numerically unstable :( Reporting the "
+                               "model to "
+                               "https://github.com/opencobra/cobrapy/issues "
+                               "will help us to fix this :)")
         LOGGER.info("found bounds infeasibility in sample, "
                     "resetting to center")
         newdir = sampler.warmup[np.random.randint(sampler.n_warmup)]
-        return _step(sampler, sampler.center, newdir - sampler.center)
+        sampler.retries += 1
+        return _step(sampler, sampler.center, newdir - sampler.center, None,
+                     tries + 1)
     return p
 
 
@@ -152,6 +162,13 @@ class HRSampler(object):
     thinning : int
         The thinning factor of the generated sampling chain. A thinning of 10
         means samples are returned every 10 steps.
+    nproj : int > 0, optional
+        How often to reproject the sampling point into the feasibility space.
+        Avoids numerical issues at the cost of lower sampling. If you observe
+        many equality constraint violations with `sampler.validate` you should
+        lower this number.
+    seed : int > 0, optional
+        The random number seed that should be used.
 
     Attributes
     ----------
@@ -162,6 +179,9 @@ class HRSampler(object):
     n_samples : int
         The total number of samples that have been generated by this
         sampler instance.
+    retries : int
+        The overall of sampling retries the sampler has observed. Larger
+        values indicate numerical instabilities.
     problem : collections.namedtuple
         A python object whose attributes define the entire sampling problem in
         matrix form. See docstring of `Problem`.
@@ -169,6 +189,8 @@ class HRSampler(object):
         A matrix of with as many columns as reactions in the model and more
         than 3 rows containing a warmup sample in each row. None if no warmup
         points have been generated yet.
+    nproj : int
+        How often to reproject the sampling point into the feasibility space.
     seed : positive integer, optional
         Sets the random number seed. Initialized to the current time stamp if
         None.
@@ -180,7 +202,7 @@ class HRSampler(object):
         the respective reverse variable.
     """
 
-    def __init__(self, model, thinning, seed=None):
+    def __init__(self, model, thinning, nproj=None, seed=None):
         """Initialize a new sampler object."""
         # This currently has to be done to reset the solver basis which is
         # required to get deterministic warmup point generation
@@ -189,7 +211,12 @@ class HRSampler(object):
             raise TypeError("sampling does not work with integer problems :(")
         self.model = model.copy()
         self.thinning = thinning
+        if nproj is None:
+            self.nproj = int(min(len(self.model.variables)**3, 1e6))
+        else:
+            self.nproj = nproj
         self.n_samples = 0
+        self.retries = 0
         self.problem = self.__build_problem()
         # Set up a map from reaction -> forward/reverse variable
         var_idx = {v: idx for idx, v in enumerate(model.variables)}
@@ -252,32 +279,54 @@ class HRSampler(object):
         if necessary).
         """
         self.n_warmup = 0
-        idx = np.hstack([self.fwd_idx, self.rev_idx])
-        self.warmup = np.zeros((len(idx), len(self.model.variables)))
+        reactions = self.model.reactions
+        self.warmup = np.zeros((2 * len(reactions), len(self.model.variables)))
         self.model.objective = Zero
-        self.model.objective.direction = "max"
-        variables = self.model.variables
-        for i in idx:
-            # Omit fixed reactions
-            if self.problem.variable_fixed[i]:
-                LOGGER.info("skipping fixed variable %s" %
-                            variables[i].name)
+        for sense in ("min", "max"):
+            self.model.objective_direction = sense
+            for i, r in enumerate(reactions):
+                variables = (self.model.variables[self.fwd_idx[i]],
+                             self.model.variables[self.rev_idx[i]])
+                # Omit fixed reactions if they are non-homogeneous
+                if r.upper_bound - r.lower_bound < bounds_tol:
+                    LOGGER.info("skipping fixed reaction %s" % r.id)
                     continue
-            self.model.objective.set_linear_coefficients({variables[i]: 1})
+                self.model.objective.set_linear_coefficients(
+                    {variables[0]: 1, variables[1]: -1})
                 self.model.slim_optimize()
                 if not self.model.solver.status == OPTIMAL:
-                LOGGER.info("can not maximize variable %s, skipping it" %
-                            variables[i].name)
+                    LOGGER.info("can not maximize reaction %s, skipping it" %
+                                r.id)
                     continue
                 primals = self.model.solver.primal_values
                 sol = [primals[v.name] for v in self.model.variables]
                 self.warmup[self.n_warmup, ] = sol
                 self.n_warmup += 1
-            # revert objective
-            self.model.objective.set_linear_coefficients({variables[i]: 0})
+                # Reset objective
+                self.model.objective.set_linear_coefficients(
+                    {variables[0]: 0, variables[1]: 0})
+        # Shrink to measure
+        self.warmup = self.warmup[0:self.n_warmup, :]
+        # Remove redundant search directions
+        keep = np.logical_not(self._is_redundant(self.warmup))
+        self.warmup = self.warmup[keep, :]
+        self.n_warmup = self.warmup.shape[0]
+
+        # Catch some special cases
+        if len(self.warmup.shape) == 1 or self.warmup.shape[0] == 1:
+            raise ValueError("Your flux cone consists only of a single point!")
+        elif self.n_warmup == 2:
+            if not self.problem.homogeneous:
+                raise ValueError("Can not sample from an inhomogenous problem"
+                                 " with only 2 search directions :(")
+            LOGGER.info("All search directions on a line, adding another one.")
+            newdir = self.warmup.T.dot([0.25, 0.25])
+            self.warmup = np.vstack([self.warmup, newdir])
+            self.n_warmup += 1
+
         # Shrink warmup points to measure
-        self.warmup = shared_np_array((self.n_warmup, len(variables)),
-                                      self.warmup[0:self.n_warmup, ])
+        self.warmup = shared_np_array(
+            (self.n_warmup, len(self.model.variables)), self.warmup)
 
     def _reproject(self, p):
         """Reproject a point into the feasibility region.
@@ -307,14 +356,28 @@ class HRSampler(object):
             new = nulls.dot(nulls.T.dot(p))
         # Projections may violate bounds
         # set to random point in space in that case
-        if any(self._bounds_dist(new) < -bounds_tol):
+        if any(new != p):
             LOGGER.info("reprojection failed in sample"
                         " %d, using random point in space" % self.n_samples)
-            idx = np.random.randint(self.n_warmup,
-                                    size=min(2, int(np.sqrt(self.n_warmup))))
-            new = self.warmup[idx, :].mean(axis=0)
+            new = self._random_point()
         return new
 
+    def _random_point(self):
+        """Find an approximately random point in the flux cone."""
+        idx = np.random.randint(self.n_warmup,
+                                size=min(2, np.ceil(np.sqrt(self.n_warmup))))
+        return self.warmup[idx, :].mean(axis=0)
+
+    def _is_redundant(self, matrix, cutoff=1.0 - feasibility_tol):
+        """Identify rdeundant rows in a matrix that can be removed."""
+        # Avoid zero variances
+        extra_col = matrix[:, 0] + 1
+        # Avoid zero rows being correlated with constant rows
+        extra_col[matrix.sum(axis=1) == 0] = 2
+        corr = np.corrcoef(np.c_[matrix, extra_col])
+        corr = np.tril(corr, -1)
+        return (np.abs(corr) > cutoff).any(axis=1)
+
     def _bounds_dist(self, p):
         """Get the lower and upper bound distances. Negative is bad."""
         prob = self.problem
@@ -445,7 +508,12 @@ class ACHRSampler(HRSampler):
     thinning : int, optional
         The thinning factor of the generated sampling chain. A thinning of 10
         means samples are returned every 10 steps.
-    seed : positive integer, optional
+    nproj : int > 0, optional
+        How often to reproject the sampling point into the feasibility space.
+        Avoids numerical issues at the cost of lower sampling. If you observe
+        many equality constraint violations with `sampler.validate` you should
+        lower this number.
+    seed : int > 0, optional
         Sets the random number seed. Initialized to the current time stamp if
         None.
 
@@ -465,9 +533,14 @@ class ACHRSampler(HRSampler):
         A matrix of with as many columns as reactions in the model and more
         than 3 rows containing a warmup sample in each row. None if no warmup
         points have been generated yet.
+    retries : int
+        The overall of sampling retries the sampler has observed. Larger
+        values indicate numerical instabilities.
     seed : positive integer, optional
         Sets the random number seed. Initialized to the current time stamp if
         None.
+    nproj : int
+        How often to reproject the sampling point into the feasibility space.
     fwd_idx : np.array
         Has one entry for each reaction in the model containing the index of
         the respective forward variable.
@@ -503,9 +576,10 @@ class ACHRSampler(HRSampler):
     .. [2] https://github.com/opencobra/cobratoolbox
     """
 
-    def __init__(self, model, thinning=100, seed=None):
+    def __init__(self, model, thinning=100, nproj=None, seed=None):
         """Initialize a new ACHRSampler."""
-        super(ACHRSampler, self).__init__(model, thinning, seed=seed)
+        super(ACHRSampler, self).__init__(model, thinning, nproj=nproj,
+                                          seed=seed)
         self.generate_fva_warmup()
         self.prev = self.center = self.warmup.mean(axis=0)
         np.random.seed(self._seed)
@@ -516,10 +590,11 @@ class ACHRSampler(HRSampler):
         delta = self.warmup[pi, ] - self.center
         self.prev = _step(self, self.prev, delta)
         if self.problem.homogeneous and (self.n_samples *
-                                         self.thinning % nproj == 0):
+                                         self.thinning % self.nproj == 0):
             self.prev = self._reproject(self.prev)
-        self.center = (self.n_samples * self.center + self.prev) / (
-                       self.n_samples + 1)
+            self.center = self._reproject(self.center)
+        self.center = ((self.n_samples * self.center) / (self.n_samples + 1) +
+                       self.prev / (self.n_samples + 1))
         self.n_samples += 1
 
     def sample(self, n, fluxes=True):
@@ -584,15 +659,17 @@ def _sample_chain(args):
         delta = sampler.warmup[pi, ] - center
 
         prev = _step(sampler, prev, delta)
-        if sampler.problem.homogeneous and (n_samples *
-                                            sampler.thinning % nproj == 0):
+        if sampler.problem.homogeneous and (
+                n_samples * sampler.thinning % sampler.nproj == 0):
             prev = sampler._reproject(prev)
+            center = sampler._reproject(center)
         if i % sampler.thinning == 0:
             samples[i//sampler.thinning - 1, ] = prev
-            center = (n_samples * center + prev) / (n_samples + 1)
+        center = ((n_samples * center) / (n_samples + 1) +
+                  prev / (n_samples + 1))
         n_samples += 1
 
-    return samples
+    return (sampler.retries, samples)
 
 
 class OptGPSampler(HRSampler):
@@ -610,7 +687,12 @@ class OptGPSampler(HRSampler):
     thinning : int, optional
         The thinning factor of the generated sampling chain. A thinning of 10
         means samples are returned every 10 steps.
-    seed : positive integer, optional
+    nproj : int > 0, optional
+        How often to reproject the sampling point into the feasibility space.
+        Avoids numerical issues at the cost of lower sampling. If you observe
+        many equality constraint violations with `sampler.validate` you should
+        lower this number.
+    seed : int > 0, optional
         Sets the random number seed. Initialized to the current time stamp if
         None.
 
@@ -630,9 +712,14 @@ class OptGPSampler(HRSampler):
         A matrix of with as many columns as reactions in the model and more
         than 3 rows containing a warmup sample in each row. None if no warmup
         points have been generated yet.
+    retries : int
+        The overall of sampling retries the sampler has observed. Larger
+        values indicate numerical instabilities.
     seed : positive integer, optional
         Sets the random number seed. Initialized to the current time stamp if
         None.
+    nproj : int
+        How often to reproject the sampling point into the feasibility space.
     fwd_idx : np.array
         Has one entry for each reaction in the model containing the index of
         the respective forward variable.
@@ -672,7 +759,8 @@ class OptGPSampler(HRSampler):
        https://doi.org/10.1371/journal.pone.0086587
     """
 
-    def __init__(self, model, processes, thinning=100, seed=None):
+    def __init__(self, model, processes, thinning=100, nproj=None,
+                 seed=None):
         """Initialize a new OptGPSampler."""
         super(OptGPSampler, self).__init__(model, thinning, seed=seed)
         self.generate_fva_warmup()
@@ -725,18 +813,19 @@ class OptGPSampler(HRSampler):
             # No with statement or starmap here since Python 2.x
             # does not support it :(
             mp = Pool(self.processes, initializer=mp_init, initargs=(self,))
-            chains = mp.map(_sample_chain, args, chunksize=1)
+            results = mp.map(_sample_chain, args, chunksize=1)
             mp.close()
             mp.join()
-            chains = np.vstack(chains)
+            chains = np.vstack([r[1] for r in results])
+            self.retries += sum(r[0] for r in results)
         else:
             mp_init(self)
-            chains = _sample_chain((n, 0))
+            results = _sample_chain((n, 0))
+            chains = results[1]
 
         # Update the global center
         self.center = (self.n_samples * self.center +
-                       n * np.atleast_2d(chains).mean(axis=0)) / (
-                       self.n_samples + n)
+                       np.atleast_2d(chains).sum(0)) / (self.n_samples + n)
         self.n_samples += n
 
         if fluxes:

cobra/flux_analysis/variability.py

@@ -2,13 +2,11 @@
 
 from __future__ import absolute_import
 
-import math
 from warnings import warn
-from itertools import chain
 
-import pandas
+from numpy import zeros
 from optlang.symbolics import Zero
-from six import iteritems
+from pandas import DataFrame
 
 from cobra.flux_analysis.loopless import loopless_fva_iter
 from cobra.flux_analysis.parsimonious import add_pfba
@@ -20,38 +18,38 @@ from cobra.util import solver as sutil
 
 def flux_variability_analysis(model, reaction_list=None, loopless=False,
                               fraction_of_optimum=1.0, pfba_factor=None):
-    """Runs flux variability analysis to find the min/max flux values for each
-    each reaction in `reaction_list`.
+    """
+    Determine the minimum and maximum possible flux value for each reaction.
 
     Parameters
     ----------
-    model : a cobra model
+    model : cobra.Model
         The model for which to run the analysis. It will *not* be modified.
     reaction_list : list of cobra.Reaction or str, optional
         The reactions for which to obtain min/max fluxes. If None will use
-        all reactions in the model.
+        all reactions in the model (default).
     loopless : boolean, optional
-        Whether to return only loopless solutions. Ignored for legacy solvers,
-        also see `Notes`.
+        Whether to return only loopless solutions. This is significantly
+        slower. Please also refer to the notes.
     fraction_of_optimum : float, optional
-        Must be <= 1.0. Requires that the objective value is at least
-        fraction * max_objective_value. A value of 0.85 for instance means that
-        the objective has to be at least at 85% percent of its maximum.
+        Must be <= 1.0. Requires that the objective value is at least the
+        fraction times maximum objective value. A value of 0.85 for instance
+        means that the objective has to be at least at 85% percent of its
+        maximum.
     pfba_factor : float, optional
-        Add additional constraint to the model that the total sum of
-        absolute fluxes must not be larger than this value times the
+        Add an additional constraint to the model that requires the total sum
+        of absolute fluxes must not be larger than this value times the
         smallest possible sum of absolute fluxes, i.e., by setting the value
-        to 1.1 then the total sum of absolute fluxes must not be more than
-        10% larger than the pfba solution. Since the pfba solution is the
-        one that optimally minimizes the total flux sum, the pfba_factor
+        to 1.1 the total sum of absolute fluxes must not be more than
+        10% larger than the pFBA solution. Since the pFBA solution is the
+        one that optimally minimizes the total flux sum, the ``pfba_factor``
         should, if set, be larger than one. Setting this value may lead to
         more realistic predictions of the effective flux bounds.
 
     Returns
     -------
     pandas.DataFrame
-        DataFrame with reaction identifier as the index columns
-
+        A data frame with reaction identifiers as the index and two columns:
         - maximum: indicating the highest possible flux
         - minimum: indicating the lowest possible flux
 
@@ -66,7 +64,7 @@ def flux_variability_analysis(model, reaction_list=None, loopless=False,
     Using the loopless option will lead to a significant increase in
     computation time (about a factor of 100 for large models). However, the
     algorithm used here (see [2]_) is still more than 1000x faster than the
-    "naive" version using `add_loopless(model)`. Also note that if you have
+    "naive" version using ``add_loopless(model)``. Also note that if you have
     included constraints that force a loop (for instance by setting all fluxes
     in a loop to be non-zero) this loop will be included in the solution.
 
@@ -85,87 +83,67 @@ def flux_variability_analysis(model, reaction_list=None, loopless=False,
     """
     if reaction_list is None:
         reaction_list = model.reactions
+    else:
+        reaction_list = model.reactions.get_by_any(reaction_list)
 
-    fva_result = _fva_optlang(model, reaction_list, fraction_of_optimum,
-                              loopless, pfba_factor)
-    return pandas.DataFrame(fva_result).T
-
-
-def _fva_optlang(model, reaction_list, fraction, loopless, pfba_factor):
-    """Helper function to perform FVA with the optlang interface.
-
-    Parameters
-    ----------
-    model : a cobra model
-    reaction_list : list of reactions
-    fraction : float, optional
-        Must be <= 1.0. Requires that the objective value is at least
-        fraction * max_objective_value. A value of 0.85 for instance means that
-        the objective has to be at least at 85% percent of its maximum.
-    loopless : boolean, optional
-        Whether to return only loopless solutions.
-    pfba_factor : float, optional
-        Add additional constraint to the model that the total sum of
-        absolute fluxes must not be larger than this value times the
-        smallest possible sum of absolute fluxes, i.e., by setting the value
-        to 1.1 then the total sum of absolute fluxes must not be more than
-        10% larger than the pfba solution. Setting this value may lead to
-        more realistic predictions of the effective flux bounds.
-
-    Returns
-    -------
-    dict
-        A dictionary containing the results.
-    """
     prob = model.problem
-    fva_results = {rxn.id: {} for rxn in reaction_list}
-    with model as m:
-        m.slim_optimize(error_value=None,
+    fva_results = DataFrame({
+        "minimum": zeros(len(reaction_list), dtype=float),
+        "maximum": zeros(len(reaction_list), dtype=float)
+    }, index=[r.id for r in reaction_list])
+    with model:
+        # Safety check before setting up FVA.
+        model.slim_optimize(error_value=None,
                             message="There is no optimal solution for the "
                                     "chosen objective!")
-        # Add objective as a variable to the model than set to zero
-        # This also uses the fraction to create the lower bound for the
-        # old objective
+        # Add the previous objective as a variable to the model then set it to
+        # zero. This also uses the fraction to create the lower/upper bound for
+        # the old objective.
+        if model.solver.objective.direction == "max":
+            fva_old_objective = prob.Variable(
+                "fva_old_objective",
+                lb=fraction_of_optimum * model.solver.objective.value)
+        else:
             fva_old_objective = prob.Variable(
-            "fva_old_objective", lb=fraction * m.solver.objective.value)
+                "fva_old_objective",
+                ub=fraction_of_optimum * model.solver.objective.value)
         fva_old_obj_constraint = prob.Constraint(
-            m.solver.objective.expression - fva_old_objective, lb=0, ub=0,
+            model.solver.objective.expression - fva_old_objective, lb=0, ub=0,
             name="fva_old_objective_constraint")
-        m.add_cons_vars([fva_old_objective, fva_old_obj_constraint])
+        model.add_cons_vars([fva_old_objective, fva_old_obj_constraint])
 
         if pfba_factor is not None:
             if pfba_factor < 1.:
-                warn('pfba_factor should be larger or equal to 1', UserWarning)
-            with m:
-                add_pfba(m, fraction_of_optimum=0)
-                ub = m.slim_optimize(error_value=None)
+                warn("The 'pfba_factor' should be larger or equal to 1.",
+                     UserWarning)
+            with model:
+                add_pfba(model, fraction_of_optimum=0)
+                ub = model.slim_optimize(error_value=None)
                 flux_sum = prob.Variable("flux_sum", ub=pfba_factor * ub)
                 flux_sum_constraint = prob.Constraint(
-                    m.solver.objective.expression - flux_sum, lb=0, ub=0,
+                    model.solver.objective.expression - flux_sum, lb=0, ub=0,
                     name="flux_sum_constraint")
-            m.add_cons_vars([flux_sum, flux_sum_constraint])
+            model.add_cons_vars([flux_sum, flux_sum_constraint])
 
-        m.objective = Zero  # This will trigger the reset as well
+        model.objective = Zero  # This will trigger the reset as well
         for what in ("minimum", "maximum"):
             sense = "min" if what == "minimum" else "max"
             for rxn in reaction_list:
-                r_id = rxn.id
-                rxn = m.reactions.get_by_id(r_id)
                 # The previous objective assignment already triggers a reset
                 # so directly update coefs here to not trigger redundant resets
                 # in the history manager which can take longer than the actual
                 # FVA for small models
-                m.solver.objective.set_linear_coefficients(
+                model.solver.objective.set_linear_coefficients(
                     {rxn.forward_variable: 1, rxn.reverse_variable: -1})
-                m.solver.objective.direction = sense
-                m.slim_optimize()
-                sutil.check_solver_status(m.solver.status)
+                model.solver.objective.direction = sense
+                model.slim_optimize()
+                sutil.check_solver_status(model.solver.status)
                 if loopless:
-                    value = loopless_fva_iter(m, rxn)
+                    value = loopless_fva_iter(model, rxn)
                 else:
-                    value = m.solver.objective.value
-                fva_results[r_id][what] = value
-                m.solver.objective.set_linear_coefficients(
+                    value = model.solver.objective.value
+                fva_results.at[rxn.id, what] = value
+                model.solver.objective.set_linear_coefficients(
                     {rxn.forward_variable: 0, rxn.reverse_variable: 0})
 
     return fva_results

cobra/test/test_flux_analysis.py

@@ -456,6 +456,13 @@ class TestCobraFluxAnalysis:
         with pytest.raises(Infeasible):
             flux_variability_analysis(infeasible_model)
 
+    def test_fva_minimization(self, model):
+        model.objective = model.reactions.EX_glc__D_e
+        model.objective_direction = 'min'
+        solution = flux_variability_analysis(model, fraction_of_optimum=.95)
+        assert solution.at['EX_glc__D_e', 'minimum'] == -10.0
+        assert solution.at['EX_glc__D_e', 'maximum'] == -9.5
+
     def test_find_blocked_reactions_solver_none(self, model):
         result = find_blocked_reactions(model, model.reactions[40:46])
         assert result == ['FRUpts2']
@@ -752,7 +759,7 @@ class TestCobraFluxSampling:
 
     def test_fixed_seed(self, model):
         s = sample(model, 1, seed=42)
-        assert numpy.allclose(s.TPI[0], 8.9516392250671544)
+        assert numpy.allclose(s.TPI[0], 9.12037487)
 
     def test_equality_constraint(self, model):
         model.reactions.ACALD.bounds = (-1.5, -1.5)
@@ -847,6 +854,48 @@ class TestCobraFluxSampling:
         proj = numpy.apply_along_axis(self.optgp._reproject, 1, s)
         assert all(self.optgp.validate(proj) == "v")
 
+    def test_complicated_model(self):
+        """Difficult model since the online mean calculation is numerically
+        unstable so many samples weakly violate the equality constraints."""
+        model = Model('flux_split')
+        reaction1 = Reaction('V1')
+        reaction2 = Reaction('V2')
+        reaction3 = Reaction('V3')
+        reaction1.lower_bound = 0
+        reaction2.lower_bound = 0
+        reaction3.lower_bound = 0
+        reaction1.upper_bound = 6
+        reaction2.upper_bound = 8
+        reaction3.upper_bound = 10
+        A = Metabolite('A')
+        reaction1.add_metabolites({A: -1})
+        reaction2.add_metabolites({A: -1})
+        reaction3.add_metabolites({A: 1})
+        model.add_reactions([reaction1])
+        model.add_reactions([reaction2])
+        model.add_reactions([reaction3])
+
+        optgp = OptGPSampler(model, 1, seed=42)
+        achr = ACHRSampler(model, seed=42)
+        optgp_samples = optgp.sample(100)
+        achr_samples = achr.sample(100)
+        assert any(optgp_samples.corr().abs() < 1.0)
+        assert any(achr_samples.corr().abs() < 1.0)
+        # > 95% are valid
+        assert(sum(optgp.validate(optgp_samples) == "v") > 95)
+        assert(sum(achr.validate(achr_samples) == "v") > 95)
+
+    def test_single_point_space(self, model):
+        """Model where constraints reduce the sampling space to one point."""
+        pfba_sol = pfba(model)
+        pfba_const = model.problem.Constraint(
+            sum(model.variables), ub=pfba_sol.objective_value)
+        model.add_cons_vars(pfba_const)
+        model.reactions.Biomass_Ecoli_core.lower_bound = \
+            pfba_sol.fluxes.Biomass_Ecoli_core
+        with pytest.raises(ValueError):
+            s = sample(model, 1)
+
 
 class TestProductionEnvelope:
     """Test the production envelope."""

cobra/test/test_solver_model.py

@@ -463,6 +463,26 @@ class TestReaction:
                 already_included_metabolite.id].expression.has(
                 -10 * reaction.reverse_variable)
 
+    def test_reaction_imul(self, model):
+        with model:
+            model.reactions.EX_glc__D_e *= 100
+            assert model.constraints.glc__D_e.expression.coeff(
+                model.variables.EX_glc__D_e) == -100
+            assert model.reactions.EX_glc__D_e.reaction == \
+                '100.0 glc__D_e <=> '
+
+        assert model.constraints.glc__D_e.expression.coeff(
+            model.variables.EX_glc__D_e) == -1
+        assert model.reactions.EX_glc__D_e.reaction == 'glc__D_e <=> '
+
+        with model:
+            model.reactions.EX_glc__D_e *= -2
+            assert model.reactions.EX_glc__D_e.bounds == (-1000.0, 10.0)
+            assert model.reactions.EX_glc__D_e.reaction == ' <=> 2.0 glc__D_e'
+
+        assert model.reactions.EX_glc__D_e.bounds == (-10, 1000.0)
+        assert model.reactions.EX_glc__D_e.reaction == 'glc__D_e <=> '
+
     # def test_pop(self, model):
     #     pgi = model.reactions.PGI
     #     g6p = model.metabolites.get_by_id("g6p_c")

cobra/util/solver.py

@@ -421,5 +421,8 @@ def assert_optimal(model, message='optimization failed'):
     message : str (optional)
         Message to for the exception if solver status was not optimal.
     """
-    if model.solver.status != optlang.interface.OPTIMAL:
-        raise OPTLANG_TO_EXCEPTIONS_DICT[model.solver.status](message)
+    status = model.solver.status
+    if status != optlang.interface.OPTIMAL:
+        exception_cls = OPTLANG_TO_EXCEPTIONS_DICT.get(
+            status, OptimizationError)
+        raise exception_cls("{} ({})".format(message, status))

debian/changelog

-python-cobra (0.11.0-1) UNRELEASED; urgency=medium
+python-cobra (0.11.3-1) UNRELEASED; urgency=medium
 
-  * New upstream version 0.11.0
+  * Team upload.
+  
+  [ Afif Elghraoui ]
+  * New upstream version
   * Update dependencies
   * Refresh patch
   * Drop obsolete patch
 
- -- Afif Elghraoui <afif@debian.org>  Sun, 04 Feb 2018 22:12:05 -0500
+  [ Andreas Tille ]
+  * Standards-Version: 4.1.4
+  * Point Vcs-fields to Salsa
+  * debhelper 11
+
+ -- Andreas Tille <tille@debian.org>  Sat, 28 Apr 2018 22:36:04 +0200
 
 python-cobra (0.5.9-1) unstable; urgency=medium
 

debian/compat

-9
+11

debian/control

 Source: python-cobra
-Section: python
-Priority: optional
 Maintainer: Debian Med Packaging Team <debian-med-packaging@lists.alioth.debian.org>
 Uploaders: Afif Elghraoui <afif@debian.org>
+Section: python
+Priority: optional
 Build-Depends:
-	debhelper (>= 9),
+	debhelper (>= 11~),
 	dh-python,
 	libglpk-dev,
 # Python2
@@ -38,11 +38,10 @@ Build-Depends:
 	python3-pytest-benchmark,
 	python-jsonschema (>> 2.5.0),
 	python3-jsonschema (>> 2.5.0),
-Standards-Version: 3.9.8
+Standards-Version: 4.1.4
+Vcs-Browser: https://salsa.debian.org/med-team/python-cobra
+Vcs-Git: https://salsa.debian.org/med-team/python-cobra.git
 Homepage: http://opencobra.github.io/cobrapy/
-Vcs-Git: https://anonscm.debian.org/git/debian-med/python-cobra.git
-Vcs-Browser: https://anonscm.debian.org/cgit/debian-med/python-cobra.git
-
 
 Package: python-cobra
 Architecture: any

documentation_builder/constraints_objectives.ipynb

 %% Cell type:markdown id: tags:
 
 # Tailored constraints, variables and objectives
 
 %% Cell type:markdown id: tags:
 
 Thanks to the use of symbolic expressions via the optlang mathematical modeling package, it is relatively straight-forward to add new variables, constraints and advanced objectives that can not easily be formulated as a combination of different reaction and their corresponding upper and lower bounds. Here we demonstrate this optlang functionality which is exposed via the `model.solver.interface`.
 
 %% Cell type:markdown id: tags:
 
 ## Constraints
 
 %% Cell type:markdown id: tags:
 
 Suppose we want to ensure that two reactions have the same flux in our model. We can add this criteria as constraint to our model using the optlang solver interface by simply defining the relevant expression as follows.
 
 %% Cell type:code id: tags:
 
 ``` python
 import cobra.test
 model = cobra.test.create_test_model('textbook')
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 same_flux = model.problem.Constraint(
     model.reactions.FBA.flux_expression - model.reactions.NH4t.flux_expression,
     lb=0,
     ub=0)
 model.add_cons_vars(same_flux)
 ```
 
 %% Cell type:markdown id: tags:
 
 The flux for our reaction of interest is obtained by the `model.reactions.FBA.flux_expression` which is simply the sum of the forward and reverse flux, i.e.,
 
 %% Cell type:code id: tags:
 
 ``` python
 model.reactions.FBA.flux_expression
 ```
 
 %% Output
 
     1.0*FBA - 1.0*FBA_reverse_84806
 
 %% Cell type:markdown id: tags:
 
 Now I can maximize growth rate whilst the fluxes of reactions 'FBA' and 'NH4t' are constrained to be (near) identical.
 
 %% Cell type:code id: tags:
 
 ``` python
 solution = model.optimize()
 print(solution.fluxes['FBA'], solution.fluxes['NH4t'],
       solution.objective_value)
 ```
 
 %% Output
 
     4.66274904774 4.66274904774 0.855110960926157
 
 %% Cell type:markdown id: tags:
 
+It is also possible to add many constraints at once. For large models, with constraints involving many reactions, the efficient way to do this is to first build a dictionary of the linear coefficients for every flux, and then add the constraint at once. For example, suppose we want to add a constrain on the sum of the absolute values of every flux in the network to be less than 100:
+
+%% Cell type:code id: tags:
+
+``` python
+coefficients = dict()
+for rxn in model.reactions:
+    coefficients[rxn.forward_variable] = 1.
+    coefficients[rxn.reverse_variable] = 1.
+constraint = model.problem.Constraint(0, lb=0, ub=100)
+model.add_cons_vars(constraint)
+model.solver.update()
+constraint.set_linear_coefficients(coefficients=coefficients)
+```
+
+%% Cell type:markdown id: tags:
+
 ## Objectives
 
 
 
 %% Cell type:markdown id: tags:
 
 Simple objective such as the maximization of the flux through one or more reactions can conveniently be done by simply
 assigning to the `model.objective` property as we have seen in previous chapters, e.g.,
 
 %% Cell type:code id: tags:
 
 ``` python
 model = cobra.test.create_test_model('textbook')
 with model:
     model.objective = {model.reactions.Biomass_Ecoli_core: 1}
     model.optimize()
     print(model.reactions.Biomass_Ecoli_core.flux)
 ```
 
 %% Output
 
     0.8739215069684307
 
 %% Cell type:markdown id: tags:
 
 The objectives mathematical expression is seen by
 
 %% Cell type:code id: tags:
 
 ``` python
 model.objective.expression
 ```
 
 %% Output
 
     -1.0*Biomass_Ecoli_core_reverse_2cdba + 1.0*Biomass_Ecoli_core
 
 %% Cell type:markdown id: tags:
 
 But suppose we need a more complicated objective, such as minimizing the Euclidean distance of the solution to the origin minus another variable, while subject to additional linear constraints. This is an objective function with both linear and quadratic components.
 
 %% Cell type:markdown id: tags:
 
 Consider the example problem:
 
 > **min** $\frac{1}{2}\left(x^2 + y^2 \right) - y$
 
 > *subject to*
 
 > $x + y = 2$
 
 > $x \ge 0$
 
 > $y \ge 0$
 
 This (admittedly very artificial) problem can be visualized graphically where the optimum is indicated by the blue dot on the line of feasible solutions.
 
 %% Cell type:code id: tags:
 
 ``` python
 %matplotlib inline
 import plot_helper
 
 plot_helper.plot_qp2()
 ```
 
 %% Output
 
 
 %% Cell type:markdown id: tags:
 
 We return to the textbook model and set the solver to one that can handle quadratic objectives such as cplex. We then add the linear constraint that the sum of our x and y reactions, that we set to FBA and NH4t, must equal 2.
 
 %% Cell type:code id: tags:
 
 ``` python
 model.solver = 'cplex'
 sum_two = model.problem.Constraint(
     model.reactions.FBA.flux_expression + model.reactions.NH4t.flux_expression,
     lb=2,
     ub=2)
 model.add_cons_vars(sum_two)
 ```
 
 %% Cell type:markdown id: tags:
 
 Next we add the quadratic objective
 
 %% Cell type:code id: tags:
 
 ``` python
 quadratic_objective = model.problem.Objective(
     0.5 * model.reactions.NH4t.flux_expression**2 + 0.5 *
     model.reactions.FBA.flux_expression**2 -
     model.reactions.FBA.flux_expression,
     direction='min')
 model.objective = quadratic_objective
 solution = model.optimize(objective_sense=None)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 print(solution.fluxes['NH4t'], solution.fluxes['FBA'])
 ```
 
 %% Output
 
     0.5 1.5
 
 %% Cell type:markdown id: tags:
 
 ## Variables
 
 %% Cell type:markdown id: tags:
 
 We can also create additional variables to facilitate studying the effects of new constraints and variables. Suppose we want to study the difference in flux between nitrogen and carbon uptake whilst we block other reactions. For this it will may help to add another variable representing this difference.
 
 %% Cell type:code id: tags:
 
 ``` python
 model = cobra.test.create_test_model('textbook')
 difference = model.problem.Variable('difference')
 ```
 
 %% Cell type:markdown id: tags:
 
 We use constraints to define what values this variable shall take
 
 %% Cell type:code id: tags:
 
 ``` python
 constraint = model.problem.Constraint(
     model.reactions.EX_glc__D_e.flux_expression -
     model.reactions.EX_nh4_e.flux_expression - difference,
     lb=0,
     ub=0)
 model.add_cons_vars([difference, constraint])
 ```
 
 %% Cell type:markdown id: tags:
 
 Now we can access that difference directly during our knock-out exploration by looking at its primal value.
 
 %% Cell type:code id: tags:
 
 ``` python
 for reaction in model.reactions[:5]:
     with model:
         reaction.knock_out()
         model.optimize()
         print(model.solver.variables.difference.primal)
 ```
 
 %% Output
 
     -5.234680806802543
     -5.2346808068025386
     -5.234680806802525
     -1.8644444444444337
     -1.8644444444444466

documentation_builder/sampling.ipynb

 %% Cell type:markdown id: tags:
 
 # Flux sampling
 
 %% Cell type:markdown id: tags:
 
 ## Basic usage
 
 %% Cell type:markdown id: tags:
 
 The easiest way to get started with flux sampling is using the `sample` function in the `flux_analysis` submodule. `sample` takes at least two arguments: a cobra model and the number of samples you want to generate.
 
 %% Cell type:code id: tags:
 
 ``` python
 from cobra.test import create_test_model
 from cobra.flux_analysis import sample
 
 model = create_test_model("textbook")
 s = sample(model, 100)
 s.head()
 ```
 
 %% Output
 
-          ACALD    ACALDt      ACKr     ACONTa     ACONTb     ACt2r       ADK1  \
-    0 -3.706944 -0.163964 -0.295823   8.975852   8.975852 -0.295823   4.847986
-    1 -1.340710 -0.175665 -0.429169  11.047827  11.047827 -0.429169   2.901598
-    2 -1.964087 -0.160334 -0.618029   9.811474   9.811474 -0.618029  17.513791
-    3 -0.838442 -0.123865 -0.376067  11.869552  11.869552 -0.376067   7.769872
-    4 -0.232088 -0.034346 -1.067684   7.972039   7.972039 -1.067684   5.114975
+          ACALD    ACALDt      ACKr     ACONTa     ACONTb     ACt2r      ADK1  \
+    0 -0.577302 -0.149662 -0.338001  10.090870  10.090870 -0.338001  0.997694
+    1 -0.639279 -0.505704 -0.031929  10.631865  10.631865 -0.031929  4.207702
+    2 -1.983410 -0.434676 -0.408318  11.046294  11.046294 -0.408318  5.510960
+    3 -1.893551 -0.618887 -0.612598   8.879426   8.879426 -0.612598  6.194372
+    4 -1.759520 -0.321021 -0.262520  10.801480  10.801480 -0.262520  4.815146
     
           AKGDH    AKGt2r    ALCD2x    ...          RPI   SUCCt2_2     SUCCt3  \
-    0  6.406533 -0.081797 -3.542980    ...    -1.649393  20.917568  20.977290
-    1  7.992916 -0.230564 -1.165045    ...    -0.066975  24.735567  24.850041
-    2  8.635576 -0.284992 -1.803753    ...    -4.075515  23.425719  23.470968
-    3  9.765178 -0.325219 -0.714577    ...    -0.838094  23.446704  23.913036
-    4  5.438125 -0.787864 -0.197742    ...    -3.109205   8.902309   9.888083
+    0  4.717467 -0.070599 -0.427639    ...    -2.255649   6.152278   6.692068
+    1  6.763224 -0.024720 -0.133575    ...    -0.769792  14.894313  14.929989
+    2  7.240802 -0.501086 -1.548735    ...    -0.088852  15.211972  15.807554
+    3  5.382222 -0.563573 -1.274664    ...    -1.728387  15.960829  17.404437
+    4  9.236588 -0.359817 -1.438499    ...    -2.840577  12.379023  13.141259
     
-            SUCDi    SUCOAS      TALA       THD2      TKT1      TKT2       TPI
-    0  744.206008 -6.406533  1.639515   1.670533  1.639515  1.635542  6.256787
-    1  710.481004 -7.992916  0.056442   9.680476  0.056442  0.052207  7.184752
-    2  696.114154 -8.635576  4.063291  52.316496  4.063291  4.058376  5.122237
-    3  595.787313 -9.765178  0.822987  36.019720  0.822987  0.816912  8.364314
-    4  584.552692 -5.438125  3.088152  12.621811  3.088152  3.079686  6.185089
+            SUCDi    SUCOAS      TALA        THD2      TKT1      TKT2       TPI
+    0  821.012351 -4.717467  2.230392  133.608893  2.230392  2.220236  5.263234
+    1  521.410118 -6.763224  0.755207   66.656770  0.755207  0.749341  7.135499
+    2  756.884622 -7.240802  0.065205   42.830676  0.065205  0.055695  8.109647
+    3  556.750972 -5.382222  1.714682   37.386748  1.714682  1.709171  7.003325
+    4  440.132011 -9.236588  2.822071    0.292885  2.822071  2.814629  6.205245
     
     [5 rows x 95 columns]
 
 %% Cell type:markdown id: tags:
 
 By default sample uses the `optgp` method based on the [method presented here](http://dx.doi.org/10.1371/journal.pone.0086587) as it is suited for larger models and can run in parallel. By default the sampler uses a single process. This can be changed by using the `processes` argument.
 
 %% Cell type:code id: tags:
 
 ``` python
 print("One process:")
 %time s = sample(model, 1000)
 print("Two processes:")
 %time s = sample(model, 1000, processes=2)
 ```
 
 %% Output
 
     One process:
-    CPU times: user 5.31 s, sys: 433 ms, total: 5.74 s
-    Wall time: 5.27 s
+    CPU times: user 7.91 s, sys: 4.09 s, total: 12 s
+    Wall time: 6.52 s
     Two processes:
-    CPU times: user 217 ms, sys: 488 ms, total: 705 ms
-    Wall time: 2.8 s
+    CPU times: user 288 ms, sys: 495 ms, total: 783 ms
+    Wall time: 3.52 s
 
 %% Cell type:markdown id: tags:
 
 Alternatively you can also user Artificial Centering Hit-and-Run for sampling by setting the method to `achr`.  `achr` does not support parallel execution but has good convergence and is almost Markovian.
 
 %% Cell type:code id: tags:
 
 ``` python
 s = sample(model, 100, method="achr")
 ```
 
 %% Cell type:markdown id: tags:
 
 In general setting up the sampler is expensive since initial search directions are generated by solving many linear programming problems. Thus, we recommend to generate as many samples as possible in one go. However, this might require finer control over the sampling procedure as described in the following section.
 
 %% Cell type:markdown id: tags:
 
 ## Advanced usage
 
 %% Cell type:markdown id: tags:
 
 ### Sampler objects
 
 %% Cell type:markdown id: tags:
 
 The sampling process can be controlled on a lower level by using the sampler classes directly.
 
 %% Cell type:code id: tags:
 
 ``` python
 from cobra.flux_analysis.sampling import OptGPSampler, ACHRSampler
 ```
 
 %% Cell type:markdown id: tags:
 
 Both sampler classes have standardized interfaces and take some additional argument. For instance the `thinning` factor. "Thinning" means only recording samples every n iterations. A higher thinning factors mean less correlated samples but also larger computation times. By default the samplers use a thinning factor of 100 which creates roughly uncorrelated samples. If you want less samples but better mixing feel free to increase this parameter. If you want to study convergence for your own model you might want to set it to 1 to obtain all iterates.
 
 %% Cell type:code id: tags:
 
 ``` python
 achr = ACHRSampler(model, thinning=10)
 ```
 
 %% Cell type:markdown id: tags:
 
 `OptGPSampler` has an additional `processes` argument specifying how many processes are used to create parallel sampling chains. This should be in the order of your CPU cores for maximum efficiency. As noted before class initialization can take up to a few minutes due to generation of initial search directions. Sampling on the other hand is quick.
 
 %% Cell type:code id: tags:
 
 ``` python
 optgp = OptGPSampler(model, processes=4)
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Sampling and validation
 
 %% Cell type:markdown id: tags:
 
 Both samplers have a sample function that generates samples from the initialized object and act like the `sample` function described above, only that this time it will only accept a single argument, the number of samples. For `OptGPSampler` the number of samples should be a multiple of the number of processes, otherwise it will be increased to the nearest multiple automatically.
 
 %% Cell type:code id: tags:
 
 ``` python
 s1 = achr.sample(100)
 
 s2 = optgp.sample(100)
 ```
 
 %% Cell type:markdown id: tags:
 
 You can call `sample` repeatedly and both samplers are optimized to generate large amount of samples without falling into "numerical traps". All sampler objects have a `validate` function in order to check if a set of points are feasible and give detailed information about feasibility violations in a form of a short code denoting feasibility. Here the short code is a combination of any of the following letters:
 
 - "v" - valid point
 - "l" - lower bound violation
 - "u" - upper bound violation
 - "e" - equality violation (meaning the point is not a steady state)
 
 For instance for a random flux distribution (should not be feasible):
 
 %% Cell type:code id: tags:
 
 ``` python
 import numpy as np
 
 bad = np.random.uniform(-1000, 1000, size=len(model.reactions))
 achr.validate(np.atleast_2d(bad))
 ```
 
 %% Output
 
-    array(['le'],
-          dtype='<U3')
+    array(['le'], dtype='<U3')
 
 %% Cell type:markdown id: tags:
 
 And for our generated samples:
 
 %% Cell type:code id: tags:
 
 ``` python
 achr.validate(s1)
 ```
 
 %% Output
 
     array(['v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v',
            'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v',
            'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v',
            'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v',
            'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v',
            'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v',
            'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v',
-           'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v'],
-          dtype='<U3')
+           'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v'], dtype='<U3')
+
+%% Cell type:markdown id: tags:
+
+Even though most models are numerically stable enought that the sampler should only generate valid samples we still urge to check this. `validate` is pretty fast and works quickly even for large models and many samples. If you find invalid samples you do not necessarily have to rerun the entire sampling but can exclude them from the sample DataFrame.
+
+%% Cell type:code id: tags:
+
+``` python
+s1_valid = s1[achr.validate(s1) == "v"]
+len(s1_valid)
+```
+
+%% Output
+
+    100
 
 %% Cell type:markdown id: tags:
 
 ### Batch sampling
 
 %% Cell type:markdown id: tags:
 
 Sampler objects are made for generating billions of samples, however using the `sample` function might quickly fill up your RAM when working with genome-scale models. Here, the `batch` method of the sampler objects might come in handy. `batch` takes two arguments, the number of samples in each batch and the number of batches. This will make sense with a small example.
 
 Let's assume we want to quantify what proportion of our samples will grow. For that we might want to generate 10 batches of 50 samples each and measure what percentage of the individual 100 samples show a growth rate larger than 0.1. Finally, we want to calculate the mean and standard deviation of those individual percentages.
 
 %% Cell type:code id: tags:
 
 ``` python
 counts = [np.mean(s.Biomass_Ecoli_core > 0.1) for s in optgp.batch(100, 10)]
 print("Usually {:.2f}% +- {:.2f}% grow...".format(
     np.mean(counts) * 100.0, np.std(counts) * 100.0))
 ```
 
 %% Output
 
-    Usually 8.70% +- 2.72% grow...
+    Usually 10.90% +- 3.83% grow...
 
 %% Cell type:markdown id: tags:
 
 ## Adding constraints
 
 %% Cell type:markdown id: tags:
 
 Flux sampling will respect additional contraints defined in the model. For instance we can add a constraint enforcing growth in asimilar manner as the section before.
 
 %% Cell type:code id: tags:
 
 ``` python
 co = model.problem.Constraint(model.reactions.Biomass_Ecoli_core.flux_expression, lb=0.1)
 model.add_cons_vars([co])
 ```
 
 %% Cell type:markdown id: tags:
 
 *Note that this is only for demonstration purposes. usually you could set the lower bound of the reaction directly instead of creating a new constraint.*
 
 %% Cell type:code id: tags:
 
 ``` python
 s = sample(model, 10)
 print(s.Biomass_Ecoli_core)
 ```
 
 %% Output
 
-    0    0.175547
-    1    0.111499
-    2    0.123073
-    3    0.151874
-    4    0.122541
-    5    0.121878
-    6    0.147333
-    7    0.106499
-    8    0.174448
-    9    0.143273
+    0    0.118106
+    1    0.120205
+    2    0.206187
+    3    0.198633
+    4    0.206575
+    5    0.119032
+    6    0.119231
+    7    0.127219
+    8    0.120086
+    9    0.182586
     Name: Biomass_Ecoli_core, dtype: float64
 
 %% Cell type:markdown id: tags:
 
 As we can see our new constraint was respected.

release-notes/0.11.1.md

+# Release notes for cobrapy 0.11.1
+
+## Fixes
+
+* Catch all `optlang` errors and re-raise them as `OptimizationError` with
+  corresponding message.

release-notes/0.11.2.md

+# Release notes for cobrapy 0.11.2
+
+## Fixes
+
+* Correctly set the constraint for a previous minimization objective in FVA.

release-notes/0.11.3.md

+# Release notes for cobrapy 0.11.3
+
+## Fixes
+
+* Improve convergence and stability of sampling.
+* Reverse inplace operations after leaving a context.
+
+## Docs
+
+* Better description of custom constraints.

release-notes/next-release.md

@@ -2,6 +2,9 @@
 
 ## Fixes
 
+* Correctly set a constraint in FVA for recording the old objective
+  value when the objective is a minimization.
+
 ## New features
 
 ## Deprecated features

setup.cfg

 [bumpversion]
-current_version = 0.11.0
+current_version = 0.11.3
 commit = True
 tag = True
 parse = (?P<major>\d+)

setup.py

@@ -36,7 +36,7 @@ except IOError:
 
 setup(
     name="cobra",
-    version="0.11.0",
+    version="0.11.3",
     packages=find_packages(),
     setup_requires=setup_requirements,
     install_requires=[
@@ -44,7 +44,7 @@ setup(
         "future",
         "swiglpk",
         "ruamel.yaml<0.15",
-        "numpy>=1.6",
+        "numpy>=1.13",
         "pandas>=0.17.0",
         "optlang>=1.2.5",
         "tabulate"