[med-svn] [python-cobra] 01/06: Imported Upstream version 0.5.6

Afif Elghraoui afif at moszumanska.debian.org
Tue Nov 29 11:47:39 UTC 2016


This is an automated email from the git hooks/post-receive script.

afif pushed a commit to branch master
in repository python-cobra.

commit 2a0a196aea6f5813f4c4d2a3728504cbd88eb91c
Author: Afif Elghraoui <afif at debian.org>
Date:   Sun Nov 27 13:01:38 2016 -0800

    Imported Upstream version 0.5.6
---
 .gitignore                                         |   3 +
 .travis.yml                                        |  96 ++-
 CONTRIBUTING.rst                                   |  44 +-
 appveyor.yml                                       |   5 +-
 cobra/VERSION                                      |   2 +-
 cobra/core/DictList.py                             |   2 +-
 cobra/core/Model.py                                |  17 +-
 cobra/core/Reaction.py                             |   3 +
 cobra/flux_analysis/parsimonious.py                |   4 +-
 cobra/io/__init__.py                               |  12 +-
 cobra/io/json.py                                   | 129 ++--
 cobra/test/__init__.py                             |  54 +-
 cobra/test/conftest.py                             |  48 ++
 cobra/test/design.py                               |  89 ---
 cobra/test/io_tests.py                             | 261 -------
 cobra/test/solvers.py                              | 292 --------
 cobra/test/test_design.py                          |  60 ++
 cobra/test/test_dictlist.py                        | 257 +++++++
 .../{flux_analysis.py => test_flux_analysis.py}    | 415 ++++++-----
 cobra/test/test_io.py                              | 213 ++++++
 .../test/{manipulation.py => test_manipulation.py} | 231 +++---
 cobra/test/test_model.py                           | 597 +++++++++++++++
 cobra/test/test_solvers.py                         | 266 +++++++
 cobra/test/unit_tests.py                           | 806 ---------------------
 config.sh                                          |  33 +-
 develop-requirements.txt                           |   1 +
 documentation_builder/conf.py                      |   2 +-
 scripts/deploy-test.sh                             |  15 -
 scripts/deploy.sh                                  |   9 +-
 tox.ini                                            |   2 +-
 30 files changed, 1990 insertions(+), 1978 deletions(-)

diff --git a/.gitignore b/.gitignore
index 05b658c..7f4e088 100644
--- a/.gitignore
+++ b/.gitignore
@@ -92,3 +92,6 @@ libglpk.a
 manylinux_builder/wheelhouse
 *~
 venv/
+.benchmarks/
+/glpk-4.60.tar.gz
+glpk-4.60
diff --git a/.travis.yml b/.travis.yml
index 345cc0a..b3a5844 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -1,43 +1,77 @@
 language: python
-sudo: false
+python: 3.5
+sudo: required
+dist: trusty
+services: docker
 cache:
   directories:
     - $HOME/.cache/pip
-python:
-  - "2.7"
-  - "3.4"
-  - "3.5"
 addons:
   apt:
     packages:
-      - gfortran
-      - libatlas-dev
-      - libatlas-base-dev
-      - liblapack-dev
-      - libgmp-dev
-      - libglpk-dev
-      - libmpfr-dev
-
-# command to install dependencies
+      - libfreetype6-dev
+      - libpng12-dev
+
 env:
-  - PIP_CACHE_DIR=$HOME/.cache/pip
+  global:
+    - secure: "hkKBaGLvoDVgktSKR3BmX+mYlGzHw9EO11MRHtiH8D9BbdygOR9p9aSV/OxkaRWhnkSP5/0SXqVgBrvU1g5OsR6cc85UQSpJ5H5jVnLoWelIbTxMCikjxDSkZlseD7ZEWrKZjRo/ZN2qym0HRWpsir3qLpl8W25xHRv/sK7Z6g8="
+    - secure: "DflyBz+QiyhlhBxn4wN00xu248EJUMjKTxUZQN6wq22qV55xO3ToGJTy9i4D6OBfZGAlSXxjjKCJ2+0sAjsghBSDEK56ud3EEg/08TIo7/T8ex/C58FsGoGFz3yDBATmquClEWN8vAMrLdxwniHmQVCBZCP/phdt5dct0AUuDc8="
+    - PLAT=x86_64
+    - UNICODE_WIDTH=32
+
+matrix:
+  exclude:
+      - python: 3.5
+  include:
+    - os: linux
+      env:
+        - MB_PYTHON_VERSION=2.7
+    - os: linux
+      env:
+        - MB_PYTHON_VERSION=3.4
+    - os: linux
+      env:
+        - MB_PYTHON_VERSION=3.5
+    - os: osx
+      language: objective-c
+      env:
+        - MB_PYTHON_VERSION=2.7
+    - os: osx
+      language: objective-c
+      env:
+        - MB_PYTHON_VERSION=3.4
+    - os: osx
+      language: objective-c
+      env:
+        - MB_PYTHON_VERSION=3.5
+
 before_install:
-  - pip install pip --upgrade
-  # These get cached
-  - pip install numpy scipy python-libsbml cython codecov jsonschema six matplotlib pandas
-  - if [[ $TRAVIS_PYTHON_VERSION == 2* ]]; then pip install lxml glpk pep8 palettable; fi
-  # Download esolver and add it to the path
-  - wget https://opencobra.github.io/pypi_cobrapy_travis/esolver.gz
-  - gzip -d esolver.gz; chmod +x esolver; export PATH=$PATH:$PWD
-  - mkdir -p $HOME/.config/matplotlib
-  - "echo 'backend: Agg' >> $HOME/.config/matplotlib/matplotlibrc"
+  - rvm get head
+  - (git clone https://github.com/matthew-brett/multibuild.git && cd multibuild && git checkout ffe5995)
+  # matplotlib non-compatible as testing runs in venv (non-framework)
+  - TEST_DEPENDS="cython codecov coverage numpy scipy python-libsbml jsonschema six pytest pytest-cov pytest-benchmark pandas tabulate"
+  - BUILD_DEPENDS="cython numpy scipy"
+  - source multibuild/common_utils.sh
+  - source multibuild/travis_steps.sh
+  - before_install
+
+before_cache:
+  - set +e
+
 install:
-  - python setup.py develop
-# # command to run tests
+  - build_wheel . $PLAT
+
 script:
-  - coverage run --source=cobra setup.py test
-  - if [[ $TRAVIS_PYTHON_VERSION == 2* ]]; then
-    pep8 cobra --exclude=oven,solvers,sbml.py --show-source;
-    fi
+  - if [[ $TRAVIS_OS_NAME == "linux" ]]; then pip install pip --upgrade; pip install pep8; pep8 cobra --exclude=oven,solvers,sbml.py --show-source; fi
+  - install_run $PLAT
+
+deploy:
+  provider: script
+  skip_cleanup: true
+  script: scripts/deploy.sh
+  on:
+    branch: master
+    tags: true
+
 after_success:
-  - codecov
+  - if [[ $TRAVIS_OS_NAME == "linux" ]]; then pip install pip --upgrade; pip install codecov; codecov; fi
diff --git a/CONTRIBUTING.rst b/CONTRIBUTING.rst
index 937460c..39e544f 100644
--- a/CONTRIBUTING.rst
+++ b/CONTRIBUTING.rst
@@ -102,8 +102,8 @@ Here's how to set up `cobrapy` for local development to contribute smaller featu
 
    Now you can make your changes locally.
 
-5. When you are done making changes, check that your changes pass pep8 and the tests with tox for the supported
-   Python versions::
+5. When you are done making changes, check that your changes pass pep8
+   and the tests with tox for the supported Python versions::
 
     (cobrapy)$ tox -e py27
     (cobrapy)$ tox -e py34
@@ -141,6 +141,46 @@ Before you submit a pull request, check that it meets these guidelines:
    Redestig. Your pull request must be approved by at least one
    reviewer before it can be merged.
 
+Unit tests and benchmarks
+-------------------------
+
+cobrapy uses `pytest <http://docs.pytest.org/en/latest/>`_ for its
+unit-tests and new features should in general always come with new
+tests that make sure that the code runs as intended. Since COBRA
+rapidly can become quite resource intensive fundamental methods such
+as model manipulation, adding and removing reactions, metabolites etc
+also must work efficiently. We use `pytest-benchmark
+<https://pytest-benchmark.readthedocs.io/en/latest/>`_ to compare
+different implementations to make sure that new code do not come with
+unacceptable increased computation time. If you add benchmarked tests,
+make sure to also include a test with and without the benchmark as we
+do not want to slow down continuous integration by running benchmarks,
+for examples, see e.g. ``test_add_metabolite`` in `test_model.py
+<cobra/test/test_model.py>`_. ``test_add_metabolite`` is the main
+test, ``test_add_metabolite_benchmark`` takes the special
+``benchmark`` fixture that enables profiling the important code
+snippet but is skipped when running::
+
+    (cobrapy)$ pytest --benchmark-skip
+
+When the test function itself is small and can safely be assumed to
+not take many resources, we can directly profile the test as in
+``test_subtract_metabolite_benchmark`` which calls
+``benchmark(self.test_subtract_metabolite, model)``.
+
+To run all tests and benchmarks do::
+
+    (cobrapy)$ pytest
+
+and to compare two implementations you may keep them in two branches
+e.g. ``old`` and ``new`` and then do::
+
+    (cobrapy)$ git checkout old
+    (cobrapy)$ pytest --benchmark-save
+    (cobrapy)$ git checkout new
+    (cobrapy)$ pytest --benchmark-compare
+
+
 Branching model
 ---------------
 
diff --git a/appveyor.yml b/appveyor.yml
index f54eb72..bb0e64e 100644
--- a/appveyor.yml
+++ b/appveyor.yml
@@ -47,12 +47,15 @@ install:
   - ps: Start-FileDownload 'https://bitbucket.org/gutworth/six/raw/default/six.py'
   - "%WITH_COMPILER% %PYTHON%/python appveyor/build_glpk.py"
   - "%PYTHON%/python -m pip install pip setuptools wheel --upgrade"
+  - "%PYTHON%/python -m pip install --upgrade pytest"
+  - "%PYTHON%/python -m pip install pytest-cov pytest-benchmark"
   - "%PYTHON%/python -m pip install Cython jsonschema twine pypandoc==1.1.3"
 
 build: off
 
 test_script:
-  - "%WITH_COMPILER% %PYTHON%/python setup.py test"
+  - "%WITH_COMPILER% %PYTHON%/python setup.py develop"
+  - "%WITH_COMPILER% %PYTHON%/python -m pytest --cov=cobra --benchmark-skip"
 
 after_test:
   - "%WITH_COMPILER% %PYTHON%/python setup.py bdist_wheel bdist_wininst"
diff --git a/cobra/VERSION b/cobra/VERSION
index 7d85683..b49b253 100644
--- a/cobra/VERSION
+++ b/cobra/VERSION
@@ -1 +1 @@
-0.5.4
+0.5.6
diff --git a/cobra/core/DictList.py b/cobra/core/DictList.py
index 2feb0ac..6d587df 100644
--- a/cobra/core/DictList.py
+++ b/cobra/core/DictList.py
@@ -79,7 +79,7 @@ class DictList(list):
                 return getattr(x, attribute)
 
         # if the search_function is a regular expression
-        if isinstance(search_function, str):
+        if isinstance(search_function, string_types):
             search_function = re.compile(search_function)
         if hasattr(search_function, "findall"):
             matches = (i for i in self
diff --git a/cobra/core/Model.py b/cobra/core/Model.py
index ab13ed0..388fb3e 100644
--- a/cobra/core/Model.py
+++ b/cobra/core/Model.py
@@ -96,17 +96,20 @@ class Model(Object):
         than deepcopy
         """
         new = self.__class__()
-        do_not_copy = {"metabolites", "reactions", "genes"}
+        do_not_copy_by_ref = {"metabolites", "reactions", "genes", "notes",
+                              "annotation"}
         for attr in self.__dict__:
-            if attr not in do_not_copy:
+            if attr not in do_not_copy_by_ref:
                 new.__dict__[attr] = self.__dict__[attr]
+        new.notes = deepcopy(self.notes)
+        new.annotation = deepcopy(self.annotation)
 
         new.metabolites = DictList()
-        do_not_copy = {"_reaction", "_model"}
+        do_not_copy_by_ref = {"_reaction", "_model"}
         for metabolite in self.metabolites:
             new_met = metabolite.__class__()
             for attr, value in iteritems(metabolite.__dict__):
-                if attr not in do_not_copy:
+                if attr not in do_not_copy_by_ref:
                     new_met.__dict__[attr] = copy(
                         value) if attr == "formula" else value
             new_met._model = new
@@ -116,18 +119,18 @@ class Model(Object):
         for gene in self.genes:
             new_gene = gene.__class__(None)
             for attr, value in iteritems(gene.__dict__):
-                if attr not in do_not_copy:
+                if attr not in do_not_copy_by_ref:
                     new_gene.__dict__[attr] = copy(
                         value) if attr == "formula" else value
             new_gene._model = new
             new.genes.append(new_gene)
 
         new.reactions = DictList()
-        do_not_copy = {"_model", "_metabolites", "_genes"}
+        do_not_copy_by_ref = {"_model", "_metabolites", "_genes"}
         for reaction in self.reactions:
             new_reaction = reaction.__class__()
             for attr, value in iteritems(reaction.__dict__):
-                if attr not in do_not_copy:
+                if attr not in do_not_copy_by_ref:
                     new_reaction.__dict__[attr] = copy(value)
             new_reaction._model = new
             new.reactions.append(new_reaction)
diff --git a/cobra/core/Reaction.py b/cobra/core/Reaction.py
index a8e5c29..d500c82 100644
--- a/cobra/core/Reaction.py
+++ b/cobra/core/Reaction.py
@@ -581,6 +581,9 @@ class Reaction(Object):
             if metabolite.charge is not None:
                 reaction_element_dict["charge"] += \
                     coefficient * metabolite.charge
+            if metabolite.elements is None:
+                raise ValueError("No elements found in metabolite %s"
+                                 % metabolite.id)
             for element, amount in iteritems(metabolite.elements):
                 reaction_element_dict[element] += coefficient * amount
         # filter out 0 values
diff --git a/cobra/flux_analysis/parsimonious.py b/cobra/flux_analysis/parsimonious.py
index 0269044..44546ea 100644
--- a/cobra/flux_analysis/parsimonious.py
+++ b/cobra/flux_analysis/parsimonious.py
@@ -47,8 +47,8 @@ def optimize_minimal_flux(cobra_model, already_irreversible=False,
                 'Minimization not supported in optimize_minimal_flux')
         optimize_kwargs.pop('objective_sense', None)
 
-    # Convert to irreversible, so all reactions will have a positive flux
-    convert_to_irreversible(cobra_model)
+    if not already_irreversible:
+        convert_to_irreversible(cobra_model)
 
     solver = solver_dict[get_solver_name() if solver is None else solver]
     lp = solver.create_problem(cobra_model, **optimize_kwargs)
diff --git a/cobra/io/__init__.py b/cobra/io/__init__.py
index d82cf88..740d4fe 100644
--- a/cobra/io/__init__.py
+++ b/cobra/io/__init__.py
@@ -6,19 +6,21 @@ from .json import load_json_model, save_json_model, to_json
 # These functions have other dependencies
 try:
     import libsbml
+    from .sbml import read_legacy_sbml
+    from .sbml import write_cobra_model_to_sbml_file as write_legacy_sbml
 except ImportError:
     warn("cobra.io.sbml requires libsbml")
     libsbml = None
-else:
-    from .sbml import read_legacy_sbml
-    from .sbml import write_cobra_model_to_sbml_file as write_legacy_sbml
+    read_legacy_sbml = None
+    write_legacy_sbml = None
 
 try:
     import scipy
+    from .mat import load_matlab_model, save_matlab_model
 except ImportError:
     warn("cobra.io.mat requires scipy")
     scipy = None
-else:
-    from .mat import load_matlab_model, save_matlab_model
+    load_matlab_model = None
+    save_matlab_model = None
 
 del libsbml, scipy, warn
diff --git a/cobra/io/json.py b/cobra/io/json.py
index fe78550..0cc8e5c 100644
--- a/cobra/io/json.py
+++ b/cobra/io/json.py
@@ -73,41 +73,47 @@ def _from_dict(obj):
     if 'reactions' not in obj:
         raise Exception('JSON object has no reactions attribute. Cannot load.')
     model = Model()
-    # add metabolites
-    new_metabolites = []
-    for metabolite in obj['metabolites']:
-        new_metabolite = Metabolite()
-        for k, v in iteritems(metabolite):
-            setattr(new_metabolite, k, v)
-        new_metabolites.append(new_metabolite)
-    model.add_metabolites(new_metabolites)
-    # add genes
-    for gene in obj['genes']:
-        new_gene = Gene(gene["id"])
-        for k, v in iteritems(gene):
-            setattr(new_gene, k, v)
-        model.genes.append(new_gene)
-    # add reactions
-    new_reactions = []
-    for reaction in obj['reactions']:
-        new_reaction = Reaction()
-        for k, v in iteritems(reaction):
-            if k == 'reversibility' or k == "reaction":
-                continue
-            elif k == 'metabolites':
-                new_reaction.add_metabolites(
-                    {model.metabolites.get_by_id(str(met)): coeff
-                     for met, coeff in iteritems(v)})
-            else:
-                setattr(new_reaction, k, v)
-        new_reactions.append(new_reaction)
-    model.add_reactions(new_reactions)
+    model.add_metabolites(
+        [metabolite_from_dict(metabolite) for metabolite in obj['metabolites']]
+    )
+    model.genes.extend([gene_from_dict(gene) for gene in obj['genes']])
+    model.add_reactions(
+        [reaction_from_dict(reaction, model) for reaction in obj['reactions']]
+    )
     for k, v in iteritems(obj):
         if k in {'id', 'name', 'notes', 'compartments', 'annotation'}:
             setattr(model, k, v)
     return model
 
 
+def reaction_from_dict(reaction, model):
+    new_reaction = Reaction()
+    for k, v in iteritems(reaction):
+        if k == 'reversibility' or k == "reaction":
+            continue
+        elif k == 'metabolites':
+            new_reaction.add_metabolites(
+                {model.metabolites.get_by_id(str(met)): coeff
+                 for met, coeff in iteritems(v)})
+        else:
+            setattr(new_reaction, k, v)
+    return new_reaction
+
+
+def gene_from_dict(gene):
+    new_gene = Gene(gene["id"])
+    for k, v in iteritems(gene):
+        setattr(new_gene, k, v)
+    return new_gene
+
+
+def metabolite_from_dict(metabolite):
+    new_metabolite = Metabolite()
+    for k, v in iteritems(metabolite):
+        setattr(new_metabolite, k, v)
+    return new_metabolite
+
+
 def _update_optional(cobra_object, new_dict, optional_attribute_dict):
     """update new_dict with optional attributes from cobra_object"""
     for key, default_value in iteritems(optional_attribute_dict):
@@ -118,42 +124,47 @@ def _update_optional(cobra_object, new_dict, optional_attribute_dict):
 
 def _to_dict(model):
     """convert the model to a dict"""
-    new_reactions = []
-    new_metabolites = []
-    new_genes = []
-    for reaction in model.reactions:
-        new_reaction = {key: _fix_type(getattr(reaction, key))
-                        for key in _REQUIRED_REACTION_ATTRIBUTES
-                        if key != "metabolites"}
-        _update_optional(reaction, new_reaction, _OPTIONAL_REACTION_ATTRIBUTES)
-        # set metabolites
-        mets = {str(met): coeff for met, coeff
-                in iteritems(reaction._metabolites)}
-        new_reaction['metabolites'] = mets
-        new_reactions.append(new_reaction)
-    for metabolite in model.metabolites:
-        new_metabolite = {key: _fix_type(getattr(metabolite, key))
-                          for key in _REQUIRED_METABOLITE_ATTRIBUTES}
-        _update_optional(metabolite, new_metabolite,
-                         _OPTIONAL_METABOLITE_ATTRIBUTES)
-        new_metabolites.append(new_metabolite)
-    for gene in model.genes:
-        new_gene = {key: str(getattr(gene, key))
-                    for key in _REQUIRED_GENE_ATTRIBUTES}
-        _update_optional(gene, new_gene, _OPTIONAL_GENE_ATTRIBUTES)
-        new_genes.append(new_gene)
-    obj = {'reactions': new_reactions,
-           'metabolites': new_metabolites,
-           'genes': new_genes,
-           'id': model.id,
-           }
-
+    obj = dict(
+        reactions=[reaction_to_dict(reaction) for reaction in model.reactions],
+        metabolites=[
+            metabolite_to_dict(metabolite) for metabolite in model.metabolites
+            ],
+        genes=[gene_to_dict(gene) for gene in model.genes],
+        id=model.id,
+    )
     _update_optional(model, obj, _OPTIONAL_MODEL_ATTRIBUTES)
     # add in the JSON version
     obj["version"] = 1
     return obj
 
 
+def gene_to_dict(gene):
+    new_gene = {key: str(getattr(gene, key))
+                for key in _REQUIRED_GENE_ATTRIBUTES}
+    _update_optional(gene, new_gene, _OPTIONAL_GENE_ATTRIBUTES)
+    return new_gene
+
+
+def metabolite_to_dict(metabolite):
+    new_metabolite = {key: _fix_type(getattr(metabolite, key))
+                      for key in _REQUIRED_METABOLITE_ATTRIBUTES}
+    _update_optional(metabolite, new_metabolite,
+                     _OPTIONAL_METABOLITE_ATTRIBUTES)
+    return new_metabolite
+
+
+def reaction_to_dict(reaction):
+    new_reaction = {key: _fix_type(getattr(reaction, key))
+                    for key in _REQUIRED_REACTION_ATTRIBUTES
+                    if key != "metabolites"}
+    _update_optional(reaction, new_reaction, _OPTIONAL_REACTION_ATTRIBUTES)
+    # set metabolites
+    mets = {str(met): coeff for met, coeff
+            in iteritems(reaction._metabolites)}
+    new_reaction['metabolites'] = mets
+    return new_reaction
+
+
 def to_json(model):
     """Save the cobra model as a json string"""
     return json.dumps(_to_dict(model), allow_nan=False)
diff --git a/cobra/test/__init__.py b/cobra/test/__init__.py
index fc59265..e30f8bb 100644
--- a/cobra/test/__init__.py
+++ b/cobra/test/__init__.py
@@ -1,31 +1,14 @@
-from __future__ import absolute_import
 from os.path import join, abspath, dirname
-import unittest as _unittest
-
+from cobra.io import read_sbml_model
+import pytest
 try:
     from cPickle import load as _load
-except:
+except ImportError:
     from pickle import load as _load
 
-from ..io import read_sbml_model
-
-
-available_tests = ['unit_tests', 'solvers', 'flux_analysis', 'io_tests',
-                   'design', 'manipulation']
-
-
 cobra_directory = abspath(join(dirname(abspath(__file__)), ".."))
 cobra_location = abspath(join(cobra_directory, ".."))
-data_directory = join(cobra_directory, "test", "data", "")
-
-salmonella_sbml = join(data_directory, "salmonella.xml")
-salmonella_pickle = join(data_directory, "salmonella.pickle")
-
-ecoli_sbml = join(data_directory, "iJO1366.xml")
-textbook_sbml = join(data_directory, "textbook.xml.gz")
-mini_sbml = join(data_directory, "mini_fbc2.xml")
-
-del abspath, join, dirname
+data_dir = join(cobra_directory, "test", "data", "")
 
 
 def create_test_model(model_name="salmonella"):
@@ -36,35 +19,24 @@ def create_test_model(model_name="salmonella"):
         path to a pickled cobra.Model
 
     """
-
     if model_name == "ecoli":
+        ecoli_sbml = join(data_dir, "iJO1366.xml")
         return read_sbml_model(ecoli_sbml)
     elif model_name == "textbook":
+        textbook_sbml = join(data_dir, "textbook.xml.gz")
         return read_sbml_model(textbook_sbml)
     elif model_name == "mini":
+        mini_sbml = join(data_dir, "mini_fbc2.xml")
         return read_sbml_model(mini_sbml)
-
-    if model_name == "salmonella":
+    elif model_name == "salmonella":
+        salmonella_pickle = join(data_dir, "salmonella.pickle")
         model_name = salmonella_pickle
     with open(model_name, "rb") as infile:
         return _load(infile)
 
 
-def create_test_suite():
-    """create a unittest.TestSuite with available tests"""
-    loader = _unittest.TestLoader()
-    suite = _unittest.TestSuite()
-    for test_name in available_tests:
-        exec("from . import " + test_name)
-        suite.addTests(loader.loadTestsFromModule(eval(test_name)))
-    return suite
-
-suite = create_test_suite()
-
-
 def test_all():
-    """###running unit tests on cobra py###"""
-    status = not _unittest.TextTestRunner(verbosity=2).run(
-        create_test_suite()
-    ).wasSuccessful()
-    return status
+    """ alias for running all unit-tests on installed cobra
+    """
+    return pytest.main(
+        ['--pyargs', 'cobra', '--benchmark-skip', '-v', '-rs']) == 0
diff --git a/cobra/test/conftest.py b/cobra/test/conftest.py
new file mode 100644
index 0000000..e9bc9b9
--- /dev/null
+++ b/cobra/test/conftest.py
@@ -0,0 +1,48 @@
+from os.path import join
+from . import create_test_model, data_dir
+import pytest
+try:
+    from cPickle import load as _load
+except ImportError:
+    from pickle import load as _load
+import json
+
+
+ at pytest.fixture(scope="session")
+def data_directory():
+    return data_dir
+
+
+ at pytest.fixture(scope="function")
+def model():
+    return create_test_model("textbook")
+
+
+ at pytest.fixture(scope="function")
+def large_model():
+    return create_test_model("ecoli")
+
+
+ at pytest.fixture(scope="function")
+def array_model():
+    return create_test_model("textbook").to_array_based_model()
+
+
+ at pytest.fixture(scope="function")
+def salmonella():
+    return create_test_model("salmonella")
+
+
+ at pytest.fixture(scope="function")
+def solved_model(data_directory):
+    model = create_test_model("textbook")
+    with open(join(data_directory, "textbook_solution.pickle"),
+              "rb") as infile:
+        model.solution = _load(infile)
+    return model
+
+
+ at pytest.fixture(scope="function")
+def fva_results(data_directory):
+    with open(join(data_directory, "textbook_fva.json"), "r") as infile:
+        return json.load(infile)
diff --git a/cobra/test/design.py b/cobra/test/design.py
deleted file mode 100644
index 0359667..0000000
--- a/cobra/test/design.py
+++ /dev/null
@@ -1,89 +0,0 @@
-from unittest import TestCase, TestLoader, TextTestRunner, skipIf
-
-import sys
-
-if __name__ == "__main__":
-    sys.path.insert(0, "../..")
-    from cobra.test import create_test_model, data_directory
-    from cobra.design import *
-    from cobra.design.design_algorithms import _add_decision_variable
-    from cobra.solvers import get_solver_name
-    sys.path.pop(0)
-else:
-    from . import create_test_model, data_directory
-    from ..design import *
-    from ..design.design_algorithms import _add_decision_variable
-    from ..solvers import get_solver_name
-
-try:
-    solver = get_solver_name(mip=True)
-except:
-    no_mip_solver = True
-else:
-    no_mip_solver = False
-
-
-class TestDesignAlgorithms(TestCase):
-    """Test functions in cobra.design"""
-
-    def test_dual(self):
-        model = create_test_model("textbook")
-        self.assertAlmostEqual(model.optimize("maximize").f, 0.874, places=3)
-        dual = dual_problem(model)
-        self.assertAlmostEqual(dual.optimize("minimize").f, 0.874, places=3)
-
-    def test_dual_integer_vars_as_lp(self):
-        model = create_test_model("textbook")
-        var = _add_decision_variable(model, "AKGDH")
-        self.assertAlmostEqual(model.optimize("maximize").f, 0.874, places=3)
-        # as lp: make integer continuous, set to 1
-        dual = dual_problem(model, "maximize", [var.id], copy=True)
-        r = dual.reactions.get_by_id(var.id)
-        r.variable_kind = "continuous"
-        r.lower_bound = r.upper_bound = 1
-        self.assertAlmostEqual(dual.optimize("minimize").f, 0.874, places=3)
-        r.lower_bound = r.upper_bound = 0
-        self.assertAlmostEqual(dual.optimize("minimize").f, 0.858, places=3)
-
-    @skipIf(no_mip_solver, "no MILP solver found")
-    def test_dual_integer_vars_as_mip(self):
-        # mip
-        model = create_test_model("textbook")
-        var = _add_decision_variable(model, "AKGDH")
-        dual = dual_problem(model, "maximize", [var.id], copy=True)
-        var_in_dual = dual.reactions.get_by_id(var.id)
-
-        # minimization, so the optimal value state is to turn off AKGDH
-        self.assertAlmostEqual(dual.optimize("minimize").f, 0.858, places=3)
-
-        # turn off AKGDH in dual
-        var_in_dual.lower_bound = var_in_dual.upper_bound = 1
-        self.assertAlmostEqual(dual.optimize("minimize").f, 0.874, places=3)
-
-        # turn on AKGDH in dual
-        var_in_dual.lower_bound = var_in_dual.upper_bound = 0
-        self.assertAlmostEqual(dual.optimize("minimize").f, 0.858, places=3)
-
-    @skipIf(no_mip_solver, "no MILP solver found")
-    def test_optknock(self):
-        model = create_test_model("textbook")
-        model.reactions.get_by_id("EX_o2_e").lower_bound = 0
-        knockable_reactions = ["ACKr", "AKGDH", "ACALD", "LDH_D"]
-        optknock_problem = set_up_optknock(model, "EX_lac__D_e",
-                                           knockable_reactions, n_knockouts=2,
-                                           copy=False)
-        solution = run_optknock(optknock_problem, tolerance_integer=1e-9)
-        self.assertIn("ACKr", solution.knockouts)
-        self.assertIn("ACALD", solution.knockouts)
-        self.assertAlmostEqual(solution.f, 17.891, places=3)
-
-# make a test suite to run all of the tests
-loader = TestLoader()
-suite = loader.loadTestsFromModule(sys.modules[__name__])
-
-
-def test_all():
-    TextTestRunner(verbosity=2).run(suite)
-
-if __name__ == "__main__":
-    test_all()
diff --git a/cobra/test/io_tests.py b/cobra/test/io_tests.py
deleted file mode 100644
index 4d7b0b5..0000000
--- a/cobra/test/io_tests.py
+++ /dev/null
@@ -1,261 +0,0 @@
-from warnings import catch_warnings, warn
-from tempfile import gettempdir
-from os import unlink
-from os.path import join, split
-from unittest import TestCase, TestLoader, TextTestRunner, skipIf, \
-    expectedFailure
-from functools import partial
-from pickle import load, dump
-from six import iteritems
-import sys
-
-if __name__ == "__main__":
-    from cobra import io
-    from cobra.test import data_directory
-else:
-    from .. import io
-    from . import data_directory
-
-libraries = ["scipy", "libsbml", "cPickle", "jsonschema"]
-for library in libraries:
-    try:
-        exec("import %s" % library)
-    except ImportError:
-        exec("%s = None" % library)
-
-with open(join(data_directory, "mini.pickle"), "rb") as infile:
-    mini_model = load(infile)
-
-
-class TestCobraIO(object):
-    def compare_models(self, model1, model2):
-        self.assertEqual(len(model1.reactions),
-                         len(model2.reactions))
-        self.assertEqual(len(model1.metabolites),
-                         len(model2.metabolites))
-        for attr in ("id", "name", "lower_bound", "upper_bound",
-                     "objective_coefficient", "gene_reaction_rule"):
-            self.assertEqual(getattr(model1.reactions[0], attr),
-                             getattr(model2.reactions[0], attr))
-            self.assertEqual(getattr(model1.reactions[10], attr),
-                             getattr(model2.reactions[10], attr))
-            self.assertEqual(getattr(model1.reactions[-1], attr),
-                             getattr(model2.reactions[-1], attr))
-        for attr in ("id", "name", "compartment", "formula", "charge"):
-            self.assertEqual(getattr(model1.metabolites[0], attr),
-                             getattr(model2.metabolites[0], attr))
-            self.assertEqual(getattr(model1.metabolites[10], attr),
-                             getattr(model2.metabolites[10], attr))
-            self.assertEqual(getattr(model1.metabolites[-1], attr),
-                             getattr(model2.metabolites[-1], attr))
-        self.assertEqual(len(model1.reactions[0].metabolites),
-                         len(model2.reactions[0].metabolites))
-        self.assertEqual(len(model1.reactions[14].metabolites),
-                         len(model2.reactions[14].metabolites))
-        self.assertEqual(len(model1.reactions[-1].metabolites),
-                         len(model2.reactions[-1].metabolites))
-        self.assertEqual(len(model1.genes), len(model2.genes))
-        # ensure they have the same solution max
-        model1.optimize()
-        model2.optimize()
-        self.assertAlmostEqual(model1.solution.f, model2.solution.f,
-                               places=3)
-        # ensure the references are correct
-        self.assertIs(model2.metabolites[0]._model, model2)
-        self.assertIs(model2.reactions[0]._model, model2)
-        self.assertIs(model2.genes[0]._model, model2)
-        self.extra_comparisons(model1, model2)
-
-    def extra_comparisons(self, model1, model2):
-        # Overriding this prevents these extra comparisons. For example, mat
-        # will need to skip this.
-        self.assertEqual(model1.compartments, model2.compartments)
-        self.assertEqual(model1.metabolites[4].annotation,
-                         model2.metabolites[4].annotation)
-        self.assertEqual(model1.reactions[4].annotation,
-                         model2.reactions[4].annotation)
-        self.assertEqual(model1.genes[4].annotation,
-                         model2.genes[4].annotation)
-        for attr in ("id", "name"):
-            self.assertEqual(getattr(model1.genes[0], attr),
-                             getattr(model2.genes[0], attr))
-            self.assertEqual(getattr(model1.genes[10], attr),
-                             getattr(model2.genes[10], attr))
-            self.assertEqual(getattr(model1.genes[-1], attr),
-                             getattr(model2.genes[-1], attr))
-        return
-
-    def test_read(self):
-        read_model = self.read_function(self.test_file)
-        self.compare_models(self.test_model, read_model)
-
-    def test_read_nonexistent(self):
-        # make sure that an error is raised when given a nonexistent file
-        self.assertRaises(IOError, self.read_function, "fake_file")
-
-    def test_write(self):
-        test_output_filename = join(gettempdir(), split(self.test_file)[-1])
-        self.write_function(self.test_model, test_output_filename)
-        reread_model = self.read_function(test_output_filename)
-        self.compare_models(self.test_model, reread_model)
-        self.validate(test_output_filename)
-        unlink(test_output_filename)
-
-    def test_write_empty(self):
-        test_output_filename = join(gettempdir(), split(self.test_file)[-1])
-        m = self.test_model.copy()
-        m.metabolites[0].charge = None
-        m.remove_reactions(list(m.reactions))
-        self.write_function(m, test_output_filename)
-        reread_model = self.read_function(test_output_filename)
-        self.assertEqual(len(reread_model.reactions), 0)
-        self.assertEqual(len(reread_model.metabolites), len(m.metabolites))
-        # ensure empty metabolite charge is read as None
-        self.assertIs(reread_model.metabolites[0].charge, None)
-        unlink(test_output_filename)
-
-    def validate(self, filename):
-        # overload if a validator exists
-        None
-
-
-class TestCobraIOSBMLfbc2(TestCase, TestCobraIO):
-    def setUp(self):
-        self.test_model = mini_model
-        self.test_file = join(data_directory, "mini_fbc2.xml")
-        self.read_function = io.read_sbml_model
-        self.write_function = io.write_sbml_model
-
-    def validate(self, filename):
-        io.sbml3.validate_sbml_model(filename)
-
-
-class TestCobraIOSBMLfbc2Gz(TestCobraIOSBMLfbc2):
-    def setUp(self):
-        self.test_model = mini_model
-        self.test_file = join(data_directory, "mini_fbc2.xml.gz")
-        self.read_function = io.read_sbml_model
-        self.write_function = io.write_sbml_model
-
-
-class TestCobraIOSBMLfbc2Bz2(TestCobraIOSBMLfbc2):
-    def setUp(self):
-        self.test_model = mini_model
-        self.test_file = join(data_directory, "mini_fbc2.xml.bz2")
-        self.read_function = io.read_sbml_model
-        self.write_function = io.write_sbml_model
-
-
-class TestCobraSBMLValidation(TestCase):
-    def test_sbml_error(self):
-        filename = join(data_directory, "invalid0.xml")
-        with self.assertRaises(io.sbml3.CobraSBMLError):
-            io.read_sbml_model(filename)
-
-    def test_bad_validation(self):
-        for i in range(3):
-            filename = join(data_directory, "invalid%d.xml" % i)
-            m, errors = io.sbml3.validate_sbml_model(filename)
-            self.assertTrue(any(len(v) >= 1 for v in iteritems(errors)))
-
-
- at skipIf(not libsbml, "libsbml required")
-class TestCobraIOSBMLfbc1(TestCase, TestCobraIO):
-    def setUp(self):
-        self.test_model = mini_model
-        self.test_file = join(data_directory, "mini_fbc1.xml")
-        self.read_function = io.read_sbml_model
-        self.write_function = partial(io.write_legacy_sbml,
-                                      use_fbc_package=True)
-
-    def extra_comparisons(self, model1, model2):
-        None
-
-    @expectedFailure
-    def test_read(self):
-        TestCobraIO.test_read(self)
-
-    @expectedFailure
-    def test_write(self):
-        TestCobraIO.test_write(self)
-
-
- at skipIf(not libsbml, "libsbml required")
-class TestCobraIOSBMLcobra(TestCase, TestCobraIO):
-    def setUp(self):
-        self.test_model = mini_model
-        self.test_file = join(data_directory, "mini_cobra.xml")
-        self.read_function = io.read_sbml_model
-        self.write_function = partial(io.write_sbml_model,
-                                      use_fbc_package=False)
-
-    def extra_comparisons(self, model1, model2):
-        None
-
-
- at skipIf(not scipy, "scipy required")
-class TestCobraIOmat(TestCase, TestCobraIO):
-    def setUp(self):
-        self.test_model = mini_model
-        self.test_file = join(data_directory, "mini.mat")
-        self.read_function = io.load_matlab_model
-        self.write_function = io.save_matlab_model
-
-    def extra_comparisons(self, model1, model2):
-        # MAT does not store gene names
-        None
-
-
-class TestCobraIOjson(TestCase, TestCobraIO):
-    def setUp(self):
-        self.test_model = mini_model
-        self.test_file = join(data_directory, "mini.json")
-        self.read_function = io.load_json_model
-        self.write_function = io.save_json_model
-
-    def validate(self, filename):
-        with open(filename, "r") as infile:
-            loaded = io.json.json.load(infile)
-        if jsonschema is None:
-            warn("jsonschema not installed")
-            return
-        jsonschema.validate(loaded, io.json.json_schema)
-
-
-class TestCobraIOPickle(TestCase, TestCobraIO):
-    def setUp(self):
-        self.test_model = mini_model
-        self.test_file = join(data_directory, "mini.pickle")
-        self.load = load
-        self.dump = dump
-
-        def read_function(filename):
-            with open(filename, "rb") as infile:
-                return self.load(infile)
-
-        def write_function(model, filename):
-            with open(filename, "wb") as outfile:
-                self.dump(model, outfile)
-
-        self.read_function = read_function
-        self.write_function = write_function
-
-
- at skipIf(not cPickle, "cPickle required")
-class TestCobraIOcPickle(TestCobraIOPickle):
-    def setUp(self):
-        TestCobraIOPickle.setUp(self)
-        self.load = cPickle.load
-        self.dump = cPickle.dump
-
-
-# make a test suite to run all of the tests
-loader = TestLoader()
-suite = loader.loadTestsFromModule(sys.modules[__name__])
-
-
-def test_all():
-    TextTestRunner(verbosity=2).run(suite)
-
-if __name__ == "__main__":
-    test_all()
diff --git a/cobra/test/solvers.py b/cobra/test/solvers.py
deleted file mode 100644
index a417e96..0000000
--- a/cobra/test/solvers.py
+++ /dev/null
@@ -1,292 +0,0 @@
-from unittest import TestCase, TestLoader, TextTestRunner, skipIf
-import sys
-# deal with absolute imports by adding the appropriate directory to the path
-if __name__ == "__main__":
-    sys.path.insert(0, "../..")
-    from cobra.test import create_test_model
-    from cobra import Model, Reaction, Metabolite
-    from cobra import solvers
-    sys.path.pop(0)  # remove the added directory to the path
-else:
-    from . import create_test_model
-    from .. import Model, Reaction, Metabolite
-    from .. import solvers
-
-try:
-    import scipy
-except:
-    scipy = None
-
-
-class TestCobraSolver(object):
-    def setUp(self):
-        self.solver = solvers.solver_dict[self.solver_name]
-        self.model = create_test_model("textbook")
-        self.old_solution = 0.8739215
-        self.infeasible_model = Model()
-        metabolite_1 = Metabolite("met1")
-        reaction_1 = Reaction("rxn1")
-        reaction_2 = Reaction("rxn2")
-        reaction_1.add_metabolites({metabolite_1: 1})
-        reaction_2.add_metabolites({metabolite_1: 1})
-        reaction_1.lower_bound = 1
-        reaction_2.upper_bound = 2
-        self.infeasible_model.add_reactions([reaction_1, reaction_2])
-
-    def test_attributes(self):
-        solver = self.solver
-        self.assertTrue(hasattr(solver, "create_problem"))
-        self.assertTrue(hasattr(solver, "solve_problem"))
-        self.assertTrue(hasattr(solver, "get_status"))
-        self.assertTrue(hasattr(solver, "get_objective_value"))
-        self.assertTrue(hasattr(solver, "format_solution"))
-        self.assertTrue(hasattr(solver, "change_variable_bounds"))
-        self.assertTrue(hasattr(solver, "change_variable_objective"))
-        self.assertTrue(hasattr(solver, "solve"))
-        self.assertTrue(hasattr(solver, "set_parameter"))
-        # self.assertTrue(hasattr(solver, "update_problem"))
-
-    def test_creation(self):
-        solver = self.solver
-        solver.create_problem(self.model)
-
-    def test_solve_feasible(self):
-        solver = self.solver
-        solution = solver.solve(self.model)
-        self.assertEqual(solution.status, "optimal")
-        self.assertAlmostEqual(self.old_solution,
-                               solution.f, places=4)
-
-    def test_solve_minimize(self):
-        solver = self.solver
-        solution = solver.solve(self.model, objective_sense='minimize')
-        self.assertEqual(solution.status, "optimal")
-        self.assertAlmostEqual(0, solution.f, places=4)
-
-    def test_low_level_control(self):
-        solver = self.solver
-        lp = solver.create_problem(self.infeasible_model)
-        solver.solve_problem(lp)
-        self.assertEqual(solver.get_status(lp), "infeasible")
-        # going to make feasible
-        solver.change_variable_bounds(lp, 0, -2., 2.)
-        solver.change_variable_bounds(lp, 1, -2., 2.)
-        solver.solve_problem(lp)
-        # should now be feasible, but obj = 0
-        self.assertEqual(solver.get_status(lp), "optimal")
-        self.assertAlmostEqual(solver.get_objective_value(lp), 0, places=4)
-        # should now have obj = 2 (maximize should be the default)
-        solver.change_variable_objective(lp, 0, 1.)
-        solver.solve_problem(lp)
-        self.assertAlmostEqual(solver.get_objective_value(lp), 2, places=4)
-        # should now solve with obj = -2
-        solver.solve_problem(lp, objective_sense="minimize")
-        self.assertAlmostEqual(solver.get_objective_value(lp), -2, places=4)
-        # should now have obj = 4
-        solver.change_variable_objective(lp, 0, 2.)
-        solver.solve_problem(lp, objective_sense="maximize")
-        self.assertAlmostEqual(solver.get_objective_value(lp), 4, places=4)
-        # make sure the solution looks good still
-        solution = solver.format_solution(lp, self.infeasible_model)
-        self.assertAlmostEqual(solution.x[0], 2, places=4)
-        self.assertAlmostEqual(solution.x[1], -2, places=4)
-        self.assertAlmostEqual(solution.x_dict["rxn1"], 2, places=4)
-        self.assertAlmostEqual(solution.x_dict["rxn2"], -2, places=4)
-
-    def test_set_objective_sense(self):
-        solver = self.solver
-        maximize = solver.create_problem(self.model,
-                                         objective_sense="maximize")
-        minimize = solver.create_problem(self.model,
-                                         objective_sense="minimize")
-        solver.solve_problem(maximize)
-        solver.solve_problem(minimize)
-        max_solution = solver.format_solution(maximize, self.model)
-        min_solution = solver.format_solution(minimize, self.model)
-        self.assertEqual(min_solution.status, "optimal")
-        self.assertAlmostEqual(0, min_solution.f, places=4)
-        self.assertAlmostEqual(self.old_solution, max_solution.f, places=4)
-        self.assertEqual(max_solution.status, "optimal")
-        # if we set minimize at creation, can we override it at solve
-        solver.solve_problem(minimize, objective_sense="maximize")
-        override_minimize = solver.format_solution(minimize, self.model)
-        self.assertAlmostEqual(max_solution.f, override_minimize.f, places=4)
-
-    def test_solve_mip(self):
-        solver = self.solver
-        if not hasattr(solver, "_SUPPORTS_MILP") or not solver._SUPPORTS_MILP:
-            self.skipTest("no milp support")
-        cobra_model = Model('MILP_implementation_test')
-        constraint = Metabolite("constraint")
-        constraint._bound = 2.5
-        x = Reaction("x")
-        x.lower_bound = 0.
-        x.objective_coefficient = 1.
-        x.add_metabolites({constraint: 2.5})
-        y = Reaction("y")
-        y.lower_bound = 0.
-        y.objective_coefficient = 1.
-        y.add_metabolites({constraint: 1.})
-        cobra_model.add_reactions([x, y])
-        float_sol = solver.solve(cobra_model)
-        # add an integer constraint
-        y.variable_kind = "integer"
-        int_sol = solver.solve(cobra_model)
-        self.assertAlmostEqual(float_sol.f, 2.5)
-        self.assertAlmostEqual(float_sol.x_dict["y"], 2.5)
-        self.assertEqual(int_sol.status, "optimal")
-        self.assertAlmostEqual(int_sol.f, 2.2)
-        self.assertAlmostEqual(int_sol.x_dict["y"], 2.0)
-
-    def test_solve_infeasible(self):
-        solver = self.solver
-        solution = solver.solve(self.infeasible_model)
-        self.assertEqual(solution.status, "infeasible")
-
-    def test_independent_creation(self):
-        solver = self.solver
-        feasible_lp = solver.create_problem(self.model)
-        infeasible_lp = solver.create_problem(self.infeasible_model)
-        solver.solve_problem(feasible_lp)
-        solver.solve_problem(infeasible_lp)
-        feasible_solution = solver.format_solution(feasible_lp, self.model)
-        infeasible_solution = solver.format_solution(infeasible_lp,
-                                                     self.infeasible_model)
-        self.assertEqual(feasible_solution.status, "optimal")
-        self.assertAlmostEqual(self.old_solution,
-                               feasible_solution.f, places=4)
-        self.assertEqual(infeasible_solution.status, "infeasible")
-
-    def test_change_coefficient(self):
-        solver = self.solver
-        c = Metabolite("c")
-        c._bound = 6
-        x = Reaction("x")
-        x.lower_bound = 1.
-        y = Reaction("y")
-        y.lower_bound = 0.
-        x.add_metabolites({c: 1})
-        z = Reaction("z")
-        z.add_metabolites({c: 1})
-        z.objective_coefficient = 1
-        m = Model("test_model")
-        m.add_reactions([x, y, z])
-        # change an existing coefficient
-        lp = solver.create_problem(m)
-        solver.solve_problem(lp)
-        sol1 = solver.format_solution(lp, m)
-        self.assertEqual(sol1.status, "optimal")
-        solver.change_coefficient(lp, 0, 0, 2)
-        solver.solve_problem(lp)
-        sol2 = solver.format_solution(lp, m)
-        self.assertEqual(sol2.status, "optimal")
-        self.assertAlmostEqual(sol1.f, 5.0)
-        self.assertAlmostEqual(sol2.f, 4.0)
-        # change a new coefficient
-        z.objective_coefficient = 0.
-        y.objective_coefficient = 1.
-        lp = solver.create_problem(m)
-        solver.change_coefficient(lp, 0, 1, 2)
-        solver.solve_problem(lp)
-        solution = solver.format_solution(lp, m)
-        self.assertEqual(solution.status, "optimal")
-        self.assertAlmostEqual(solution.x_dict["y"], 2.5)
-
-    def test_inequality(self):
-        # The space enclosed by the constraints is a 2D triangle with
-        # vertexes as (3, 0), (1, 2), and (0, 1)
-        solver = self.solver
-        # c1 encodes y - x > 1 ==> y > x - 1
-        # c2 encodes y + x < 3 ==> y < 3 - x
-        c1 = Metabolite("c1")
-        c2 = Metabolite("c2")
-        x = Reaction("x")
-        x.lower_bound = 0
-        y = Reaction("y")
-        y.lower_bound = 0
-        x.add_metabolites({c1: -1, c2: 1})
-        y.add_metabolites({c1: 1, c2: 1})
-        c1._bound = 1
-        c1._constraint_sense = "G"
-        c2._bound = 3
-        c2._constraint_sense = "L"
-        m = Model()
-        m.add_reactions([x, y])
-        # test that optimal values are at the vertices
-        m.objective = "x"
-        self.assertAlmostEqual(solver.solve(m).f, 1.0)
-        self.assertAlmostEqual(solver.solve(m).x_dict["y"], 2.0)
-        m.objective = "y"
-        self.assertAlmostEqual(solver.solve(m).f, 3.0)
-        self.assertAlmostEqual(solver.solve(m, objective_sense="minimize").f,
-                               1.0)
-
-    @skipIf(scipy is None, "scipy required for quadratic objectives")
-    def test_quadratic(self):
-        solver = self.solver
-        if not hasattr(solver, "set_quadratic_objective"):
-            self.skipTest("no qp support")
-        c = Metabolite("c")
-        c._bound = 2
-        x = Reaction("x")
-        x.objective_coefficient = -0.5
-        x.lower_bound = 0.
-        y = Reaction("y")
-        y.objective_coefficient = -0.5
-        y.lower_bound = 0.
-        x.add_metabolites({c: 1})
-        y.add_metabolites({c: 1})
-        m = Model()
-        m.add_reactions([x, y])
-        lp = self.solver.create_problem(m)
-        quadratic_obj = scipy.sparse.eye(2) * 2
-        solver.set_quadratic_objective(lp, quadratic_obj)
-        solver.solve_problem(lp, objective_sense="minimize")
-        solution = solver.format_solution(lp, m)
-        self.assertEqual(solution.status, "optimal")
-        # Respecting linear objectives also makes the objective value 1.
-        self.assertAlmostEqual(solution.f, 1.)
-        self.assertAlmostEqual(solution.x_dict["y"], 1.)
-        self.assertAlmostEqual(solution.x_dict["y"], 1.)
-        # When the linear objectives are removed the objective value is 2.
-        solver.change_variable_objective(lp, 0, 0.)
-        solver.change_variable_objective(lp, 1, 0.)
-        solver.solve_problem(lp, objective_sense="minimize")
-        solution = solver.format_solution(lp, m)
-        self.assertEqual(solution.status, "optimal")
-        self.assertAlmostEqual(solution.f, 2.)
-        # test quadratic from solve function
-        solution = solver.solve(m, quadratic_component=quadratic_obj,
-                                objective_sense="minimize")
-        self.assertEqual(solution.status, "optimal")
-        self.assertAlmostEqual(solution.f, 1.)
-        c._bound = 6
-        z = Reaction("z")
-        x.objective_coefficient = 0.
-        y.objective_coefficient = 0.
-        z.lower_bound = 0.
-        z.add_metabolites({c: 1})
-        m.add_reaction(z)
-        solution = solver.solve(m, quadratic_component=scipy.sparse.eye(3),
-                                objective_sense="minimize")
-        # should be 12 not 24 because 1/2 (V^T Q V)
-        self.assertEqual(solution.status, "optimal")
-        self.assertAlmostEqual(solution.f, 6)
-        self.assertAlmostEqual(solution.x_dict["x"], 2, places=6)
-        self.assertAlmostEqual(solution.x_dict["y"], 2, places=6)
-        self.assertAlmostEqual(solution.x_dict["z"], 2, places=6)
-
-for solver_name in solvers.solver_dict:
-    exec('class %sTester(TestCobraSolver, TestCase): None' % solver_name)
-    exec('%sTester.solver_name = "%s"' % (solver_name, solver_name))
-
-# make a test suite to run all of the tests
-loader = TestLoader()
-suite = loader.loadTestsFromModule(sys.modules[__name__])
-
-
-def test_all():
-    TextTestRunner(verbosity=2).run(suite)
-
-if __name__ == "__main__":
-    test_all()
diff --git a/cobra/test/test_design.py b/cobra/test/test_design.py
new file mode 100644
index 0000000..4c65565
--- /dev/null
+++ b/cobra/test/test_design.py
@@ -0,0 +1,60 @@
+import pytest
+
+from cobra.design import *
+from cobra.design.design_algorithms import _add_decision_variable
+from cobra.solvers import get_solver_name
+from .conftest import model
+
+try:
+    solver = get_solver_name(mip=True)
+except ImportError:
+    no_mip_solver = True
+else:
+    no_mip_solver = False
+
+
+class TestDesignAlgorithms:
+    """Test functions in cobra.design"""
+    def test_dual(self, model):
+        assert abs(model.optimize("maximize").f - 0.874) < 0.001
+        dual = dual_problem(model)
+        assert abs(dual.optimize("minimize").f - 0.874) < 0.001
+
+    def test_dual_integer_vars_as_lp(self, model):
+        var = _add_decision_variable(model, "AKGDH")
+        assert abs(model.optimize("maximize").f - 0.874) < 0.001
+        # as lp: make integer continuous, set to 1
+        dual = dual_problem(model, "maximize", [var.id], copy=True)
+        r = dual.reactions.get_by_id(var.id)
+        r.variable_kind = "continuous"
+        r.lower_bound = r.upper_bound = 1
+        assert abs(dual.optimize("minimize").f - 0.874) < 0.001
+        r.lower_bound = r.upper_bound = 0
+        assert abs(dual.optimize("minimize").f - 0.858) < 0.001
+
+    @pytest.mark.skipif(no_mip_solver, reason="no MILP solver found")
+    def test_dual_integer_vars_as_mip(self, model):
+        # mip
+        var = _add_decision_variable(model, "AKGDH")
+        dual = dual_problem(model, "maximize", [var.id], copy=True)
+        var_in_dual = dual.reactions.get_by_id(var.id)
+        # minimization, so the optimal value state is to turn off AKGDH
+        assert abs(dual.optimize("minimize").f - 0.858) < 0.001
+        # turn off AKGDH in dual
+        var_in_dual.lower_bound = var_in_dual.upper_bound = 1
+        assert abs(dual.optimize("minimize").f - 0.874) < 0.001
+        # turn on AKGDH in dual
+        var_in_dual.lower_bound = var_in_dual.upper_bound = 0
+        assert abs(dual.optimize("minimize").f - 0.858) < 0.001
+
+    @pytest.mark.skipif(no_mip_solver, reason="no MILP solver found")
+    def test_optknock(self, model):
+        model.reactions.get_by_id("EX_o2_e").lower_bound = 0
+        knockable_reactions = ["ACKr", "AKGDH", "ACALD", "LDH_D"]
+        optknock_problem = set_up_optknock(model, "EX_lac__D_e",
+                                           knockable_reactions, n_knockouts=2,
+                                           copy=False)
+        solution = run_optknock(optknock_problem, tolerance_integer=1e-9)
+        assert "ACKr" in solution.knockouts
+        assert "ACALD" in solution.knockouts
+        assert abs(solution.f - 17.891) < 0.001
diff --git a/cobra/test/test_dictlist.py b/cobra/test/test_dictlist.py
new file mode 100644
index 0000000..e7e072b
--- /dev/null
+++ b/cobra/test/test_dictlist.py
@@ -0,0 +1,257 @@
+import pytest
+from copy import deepcopy, copy
+import re
+from cobra import DictList, Object
+from pickle import loads, dumps, HIGHEST_PROTOCOL
+
+
+ at pytest.fixture(scope="function")
+def dict_list():
+    obj = Object("test1")
+    test_list = DictList()
+    test_list.append(obj)
+    return obj, test_list
+
+
+class TestDictList:
+    def test_contains(self, dict_list):
+        obj, test_list = dict_list
+        assert obj in test_list
+        assert obj.id in test_list
+        assert Object("not_in") not in test_list
+        assert 'not_in' not in test_list
+
+    def test_index(self, dict_list):
+        obj, test_list = dict_list
+        assert test_list.index("test1") == 0
+        assert test_list.index(obj) == 0
+        with pytest.raises(ValueError):
+            test_list.index("f")
+        with pytest.raises(ValueError):
+            test_list.index(Object("f"))
+        # ensure trying to index with an object that is a different object
+        # also raises an error
+        with pytest.raises(ValueError):
+            test_list.index(Object("test1"))
+
+    def test_independent(self):
+        a = DictList([Object("o1"), Object("o2")])
+        b = DictList()
+        assert "o1" in a
+        assert "o1" not in b
+        b.append(Object("o3"))
+        assert "o3" not in a
+        assert "o3" in b
+
+    def test_append(self, dict_list):
+        obj, test_list = dict_list
+        obj2 = Object("test2")
+        test_list.append(obj2)
+        with pytest.raises(ValueError):
+            test_list.append(Object("test1"))
+        assert test_list.index(obj2) == 1
+        assert test_list[1] == obj2
+        assert test_list.get_by_id("test2") is obj2
+        assert len(test_list) == 2
+
+    def test_insert(self, dict_list):
+        obj, test_list = dict_list
+        obj2 = Object("a")
+        test_list.insert(0, obj2)
+        assert test_list.index(obj2) == 0
+        assert test_list.index("test1") == 1
+        assert test_list.get_by_id("a") is obj2
+        assert len(test_list) == 2
+        with pytest.raises(ValueError):
+            test_list.append(obj2)
+
+    def test_extend(self, dict_list):
+        obj, test_list = dict_list
+        obj_list = [Object("test%d" % (i)) for i in range(2, 10)]
+        test_list.extend(obj_list)
+        assert test_list[1].id == "test2"
+        assert test_list.get_by_id("test2") == obj_list[0]
+        assert test_list[8].id == "test9"
+        assert len(test_list) == 9
+        with pytest.raises(ValueError):
+            test_list.extend([Object("test1")])
+        # Even if the object is unique, if it is present twice in the new
+        # list, it should still raise an exception
+        with pytest.raises(ValueError):
+            test_list.extend([Object("testd"), Object("testd")])
+
+    def test_iadd(self, dict_list):
+        obj, test_list = dict_list
+        obj_list = [Object("test%d" % (i)) for i in range(2, 10)]
+        test_list += obj_list
+        assert test_list[1].id == "test2"
+        assert test_list.get_by_id("test2") == obj_list[0]
+        assert test_list[8].id == "test9"
+        assert len(test_list) == 9
+
+    def test_add(self, dict_list):
+        obj, test_list = dict_list
+        obj_list = [Object("test%d" % (i)) for i in range(2, 10)]
+        sum = test_list + obj_list
+        assert sum is not test_list
+        assert sum is not obj_list
+        assert test_list[0].id == "test1"
+        assert sum[1].id == "test2"
+        assert sum.get_by_id("test2") == obj_list[0]
+        assert sum[8].id == "test9"
+        assert len(test_list) == 1
+        assert len(sum) == 9
+
+    def test_init_copy(self, dict_list):
+        obj, test_list = dict_list
+        test_list.append(Object("test2"))
+        copied = DictList(test_list)
+        assert test_list is not copied
+        assert isinstance(copied, test_list.__class__)
+        assert len(test_list) == len(copied)
+        for i, v in enumerate(test_list):
+            assert test_list[i].id == copied[i].id
+            assert i == copied.index(v.id)
+            assert test_list[i] is copied[i]
+            assert v is copied.get_by_id(v.id)
+
+    def test_slice(self, dict_list):
+        obj, test_list = dict_list
+        test_list.append(Object("test2"))
+        test_list.append(Object("test3"))
+        sliced = test_list[:-1]
+        assert test_list is not sliced
+        assert isinstance(sliced, test_list.__class__)
+        assert len(test_list) == len(sliced) + 1
+        for i, v in enumerate(sliced):
+            assert test_list[i].id == sliced[i].id
+            assert i == sliced.index(v.id)
+            assert test_list[i] is sliced[i]
+            assert test_list[i] is sliced.get_by_id(v.id)
+
+    def test_copy(self, dict_list):
+        obj, test_list = dict_list
+        test_list.append(Object("test2"))
+        copied = copy(test_list)
+        assert test_list is not copied
+        assert isinstance(copied, test_list.__class__)
+        assert len(test_list) == len(copied)
+        for i, v in enumerate(test_list):
+            assert test_list[i].id == copied[i].id
+            assert i == copied.index(v.id)
+            assert test_list[i] is copied[i]
+            assert v is copied.get_by_id(v.id)
+
+    def test_deepcopy(self, dict_list):
+        obj, test_list = dict_list
+        test_list.append(Object("test2"))
+        copied = deepcopy(test_list)
+        assert test_list is not copied
+        assert isinstance(copied, test_list.__class__)
+        assert len(test_list) == len(copied)
+        for i, v in enumerate(test_list):
+            assert test_list[i].id == copied[i].id
+            assert i == copied.index(v.id)
+            assert test_list[i] is not copied[i]
+            assert v is not copied.get_by_id(v.id)
+
+    def test_pickle(self, dict_list):
+        obj, test_list = dict_list
+        test_list.append(Object("test2"))
+        for protocol in range(HIGHEST_PROTOCOL):
+            pickle_str = dumps(test_list, protocol=protocol)
+            copied = loads(pickle_str)
+            assert test_list is not copied
+            assert isinstance(copied, test_list.__class__)
+            assert len(test_list) == len(copied)
+            for i, v in enumerate(test_list):
+                assert test_list[i].id == copied[i].id
+                assert i == copied.index(v.id)
+                assert test_list[i] is not copied[i]
+                assert v is not copied.get_by_id(v.id)
+
+    def test_query(self, dict_list):
+        obj, test_list = dict_list
+        obj2 = Object("test2")
+        obj2.name = "foobar1"
+        test_list.append(obj2)
+        result = test_list.query("test1")  # matches only test1
+        assert len(result) == 1
+        result = test_list.query(u"test1")  # matches with unicode
+        assert len(result) == 1
+        assert result[0] == obj
+        result = test_list.query("foo", "name")  # matches only test2
+        assert len(result) == 1
+        assert result[0] == obj2
+        result = test_list.query("test")  # matches test1 and test2
+        assert len(result) == 2
+        # test with a regular expression
+        result = test_list.query(re.compile("test[0-9]"))
+        assert len(result) == 2
+        result = test_list.query(re.compile("test[29]"))
+        assert len(result) == 1
+        # test query of name
+        result = test_list.query(re.compile("foobar."), "name")
+        assert len(result) == 1
+
+    def test_removal(self):
+        obj_list = DictList(Object("test%d" % (i)) for i in range(2, 10))
+        del obj_list[3]
+        assert "test5" not in obj_list
+        assert obj_list.index(obj_list[-1]) == len(obj_list) - 1
+        assert len(obj_list) == 7
+        del obj_list[3:5]
+        assert "test6" not in obj_list
+        assert "test7" not in obj_list
+        assert obj_list.index(obj_list[-1]) == len(obj_list) - 1
+        assert len(obj_list) == 5
+        removed = obj_list.pop(1)
+        assert obj_list.index(obj_list[-1]) == len(obj_list) - 1
+        assert removed.id == "test3"
+        assert "test3" not in obj_list
+        assert len(obj_list) == 4
+        removed = obj_list.pop()
+        assert removed.id == "test9"
+        assert removed.id not in obj_list
+        assert len(obj_list) == 3
+
+    def test_set(self):
+        obj_list = DictList(Object("test%d" % (i)) for i in range(10))
+        obj_list[4] = Object("testa")
+        assert obj_list.index("testa") == 4
+        assert obj_list[4].id == "testa"
+        obj_list[5:7] = [Object("testb"), Object("testc")]
+        assert obj_list.index("testb") == 5
+        assert obj_list[5].id == "testb"
+        assert obj_list.index("testc") == 6
+        assert obj_list[6].id == "testc"
+        # Even if the object is unique, if it is present twice in the new
+        # list, it should still raise an exception
+        with pytest.raises(ValueError):
+            obj_list.__setitem__(slice(5, 7),
+                                 [Object("testd"), Object("testd")])
+
+    def test_sort_and_reverse(self):
+        dl = DictList(Object("test%d" % (i)) for i in reversed(range(10)))
+        assert dl[0].id == "test9"
+        dl.sort()
+        assert len(dl) == 10
+        assert dl[0].id == "test0"
+        assert dl.index("test0") == 0
+        dl.reverse()
+        assert dl[0].id == "test9"
+        assert dl.index("test0") == 9
+
+    def test_dir(self, dict_list):
+        obj, test_list = dict_list
+        """makes sure tab complete will work"""
+        attrs = dir(test_list)
+        assert "test1" in attrs
+        assert "_dict" in attrs  # attribute of DictList
+
+    def test_union(self, dict_list):
+        obj, test_list = dict_list
+        test_list.union([Object("test1"), Object("test2")])
+        # should only add 1 element
+        assert len(test_list) == 2
+        assert test_list.index("test2") == 1
diff --git a/cobra/test/flux_analysis.py b/cobra/test/test_flux_analysis.py
similarity index 50%
rename from cobra/test/flux_analysis.py
rename to cobra/test/test_flux_analysis.py
index 8c081a2..36b6ca9 100644
--- a/cobra/test/flux_analysis.py
+++ b/cobra/test/test_flux_analysis.py
@@ -1,48 +1,32 @@
-from unittest import TestCase, TestLoader, TextTestRunner, skipIf
-
-from warnings import warn
+import pytest
 import sys
-from os.path import join
 from os import name
-from json import load
 from contextlib import contextmanager
-import pickle
 import re
-
 from six import iteritems, StringIO
+from cobra.core import Model, Reaction, Metabolite
+from cobra.solvers import solver_dict, get_solver_name
+from cobra.flux_analysis import *
+from cobra.solvers import SolverNotFound
+from .conftest import model, large_model, solved_model, fva_results
 
 try:
     import numpy
-except:
+except ImportError:
     numpy = None
 try:
     import matplotlib
-except:
+except ImportError:
     matplotlib = None
 try:
     import pandas
-except:
+except ImportError:
     pandas = None
 try:
     import tabulate
-except:
+except ImportError:
     tabulate = None
 
-if __name__ == "__main__":
-    sys.path.insert(0, "../..")
-    from cobra.test import create_test_model, data_directory
-    from cobra import Model, Reaction, Metabolite
-    from cobra.manipulation import initialize_growth_medium
-    from cobra.solvers import solver_dict, get_solver_name
-    from cobra.flux_analysis import *
-    sys.path.pop(0)
-else:
-    from . import create_test_model, data_directory
-    from .. import Model, Reaction, Metabolite
-    from ..manipulation import initialize_growth_medium
-    from ..solvers import solver_dict, get_solver_name
-    from ..flux_analysis import *
-
 
 @contextmanager
 def captured_output():
@@ -56,124 +40,144 @@ def captured_output():
         sys.stdout, sys.stderr = old_out, old_err
 
 
-class TestCobraFluxAnalysis(TestCase):
+class TestCobraFluxAnalysis:
     """Test the simulation functions in cobra.flux_analysis"""
 
-    def setUp(self):
-        pass
-
-    def test_pFBA(self):
-        model = create_test_model("textbook")
-        for solver in solver_dict:
+    @pytest.mark.parametrize("solver", list(solver_dict))
+    def test_pfba_benchmark(self, large_model, benchmark, solver):
+        benchmark(optimize_minimal_flux, large_model, solver=solver)
+
+    @pytest.mark.parametrize("solver", list(solver_dict))
+    def test_pfba(self, model, solver):
+        optimize_minimal_flux(model, solver=solver)
+        abs_x = [abs(i) for i in model.solution.x]
+        assert model.solution.status == "optimal"
+        assert abs(model.solution.f - 0.8739) < 0.001
+        assert abs(sum(abs_x) - 518.4221) < 0.001
+
+        # Test desired_objective_value
+        desired_objective = 0.8
+        optimize_minimal_flux(model, solver=solver,
+                              desired_objective_value=desired_objective)
+        abs_x = [abs(i) for i in model.solution.x]
+        assert model.solution.status == "optimal"
+        assert abs(model.solution.f - desired_objective) < 0.001
+        assert abs(sum(abs_x) - 476.1594) < 0.001
+
+        # Test fraction_of_optimum
+        optimize_minimal_flux(model, solver=solver,
+                              fraction_of_optimum=0.95)
+        abs_x = [abs(i) for i in model.solution.x]
+        assert model.solution.status == "optimal"
+        assert abs(model.solution.f - 0.95 * 0.8739) < 0.001
+        assert abs(sum(abs_x) - 493.4400) < 0.001
+
+        # Make sure the model works for non-unity objective values
+        model.reactions.Biomass_Ecoli_core.objective_coefficient = 2
+        optimize_minimal_flux(model, solver=solver)
+        assert abs(model.solution.f - 2 * 0.8739) < 0.001
+        model.reactions.Biomass_Ecoli_core.objective_coefficient = 1
+
+        # Test some erroneous inputs -- multiple objectives
+        model.reactions.ATPM.objective_coefficient = 1
+        with pytest.raises(ValueError):
             optimize_minimal_flux(model, solver=solver)
-            abs_x = [abs(i) for i in model.solution.x]
-            self.assertEqual(model.solution.status, "optimal")
-            self.assertAlmostEqual(model.solution.f, 0.8739, places=3)
-            self.assertAlmostEqual(sum(abs_x), 518.4221, places=3)
+        model.reactions.ATPM.objective_coefficient = 0
 
-            # Test desired_objective_value
-            desired_objective = 0.8
+        # Minimization of objective
+        with pytest.raises(ValueError):
             optimize_minimal_flux(model, solver=solver,
-                                  desired_objective_value=desired_objective)
-            abs_x = [abs(i) for i in model.solution.x]
-            self.assertEqual(model.solution.status, "optimal")
-            self.assertAlmostEqual(model.solution.f, desired_objective,
-                                   places=3)
-            self.assertAlmostEqual(sum(abs_x), 476.1594, places=3)
-
-            # Test fraction_of_optimum
-            optimize_minimal_flux(model, solver=solver,
-                                  fraction_of_optimum=0.95)
-            abs_x = [abs(i) for i in model.solution.x]
-            self.assertEqual(model.solution.status, "optimal")
-            self.assertAlmostEqual(model.solution.f, 0.95*0.8739, places=3)
-            self.assertAlmostEqual(sum(abs_x), 493.4400, places=3)
-
-            # Make sure the model works for non-unity objective values
-            model.reactions.Biomass_Ecoli_core.objective_coefficient = 2
+                                  objective_sense='minimize')
+
+        # Infeasible solution
+        atpm = float(model.reactions.ATPM.lower_bound)
+        model.reactions.ATPM.lower_bound = 500
+        with pytest.raises(ValueError):
             optimize_minimal_flux(model, solver=solver)
-            self.assertAlmostEqual(model.solution.f, 2*0.8739, places=3)
-            model.reactions.Biomass_Ecoli_core.objective_coefficient = 1
-
-            # Test some erroneous inputs -- multiple objectives
-            model.reactions.ATPM.objective_coefficient = 1
-            with self.assertRaises(ValueError):
-                optimize_minimal_flux(model, solver=solver)
-            model.reactions.ATPM.objective_coefficient = 0
-
-            # Minimization of objective
-            with self.assertRaises(ValueError):
-                optimize_minimal_flux(model, solver=solver,
-                                      objective_sense='minimize')
-
-            # Infeasible solution
-            atpm = float(model.reactions.ATPM.lower_bound)
-            model.reactions.ATPM.lower_bound = 500
-            with self.assertRaises(ValueError):
-                optimize_minimal_flux(model, solver=solver)
-            model.reactions.ATPM.lower_bound = atpm
-
-    def test_single_gene_deletion_fba(self):
-        cobra_model = create_test_model("textbook")
+        model.reactions.ATPM.lower_bound = atpm
+
+    def test_single_gene_deletion_fba_benchmark(self, large_model, benchmark):
+        genes = ['b0511', 'b2521', 'b0651', 'b2502', 'b3132', 'b1486', 'b3384',
+                 'b4321', 'b3428', 'b2789', 'b0052', 'b0115',
+                 'b2167', 'b0759', 'b3389', 'b4031', 'b3916', 'b2374', 'b0677',
+                 'b2202']
+        benchmark(single_gene_deletion, large_model, gene_list=genes)
+
+    def test_single_gene_deletion_fba(self, model):
         # expected knockouts for textbook model
         growth_dict = {"b0008": 0.87, "b0114": 0.80, "b0116": 0.78,
                        "b2276": 0.21, "b1779": 0.00}
-
-        rates, statuses = single_gene_deletion(cobra_model,
+        rates, statuses = single_gene_deletion(model,
                                                gene_list=growth_dict.keys(),
                                                method="fba")
         for gene, expected_value in iteritems(growth_dict):
-            self.assertEqual(statuses[gene], 'optimal')
-            self.assertAlmostEqual(rates[gene], expected_value, places=2)
+            assert statuses[gene] == 'optimal'
+            assert abs(rates[gene] - expected_value) < 0.01
+
+    def test_single_gene_deletion_moma_benchmark(self, large_model, benchmark):
+        try:
+            get_solver_name(qp=True)
+        except SolverNotFound:
+            pytest.skip("no qp support")
+        genes = ['b1764', 'b0463', 'b1779', 'b0417']
+        benchmark(single_gene_deletion, large_model, gene_list=genes,
+                  method="moma")
 
-    def test_single_gene_deletion_moma(self):
-        # MOMA requires a QP solver
+    def test_single_gene_deletion_moma(self, model):
         try:
             get_solver_name(qp=True)
-        except:
-            self.skipTest("no qp support")
+        except SolverNotFound:
+            pytest.skip("no qp support")
 
-        cobra_model = create_test_model("textbook")
         # expected knockouts for textbook model
         growth_dict = {"b0008": 0.87, "b0114": 0.71, "b0116": 0.56,
                        "b2276": 0.11, "b1779": 0.00}
 
-        rates, statuses = single_gene_deletion(cobra_model,
+        rates, statuses = single_gene_deletion(model,
                                                gene_list=growth_dict.keys(),
                                                method="moma")
         for gene, expected_value in iteritems(growth_dict):
-            self.assertEqual(statuses[gene], 'optimal')
-            self.assertAlmostEqual(rates[gene], expected_value, places=2)
+            assert statuses[gene] == 'optimal'
+            assert abs(rates[gene] - expected_value) < 0.01
 
-    def test_single_reaction_deletion(self):
-        cobra_model = create_test_model("textbook")
+    def test_single_gene_deletion_benchmark(self, large_model, benchmark):
+        reactions = ['CDPMEK', 'PRATPP', 'HISTD', 'PPCDC']
+        benchmark(single_reaction_deletion, large_model,
+                  reaction_list=reactions)
+
+    def test_single_reaction_deletion(self, model):
         expected_results = {'FBA': 0.70404, 'FBP': 0.87392, 'CS': 0,
                             'FUM': 0.81430, 'GAPD': 0, 'GLUDy': 0.85139}
 
         results, status = single_reaction_deletion(
-            cobra_model, reaction_list=expected_results.keys())
-        self.assertEqual(len(results), 6)
-        self.assertEqual(len(status), 6)
+            model, reaction_list=expected_results.keys())
+        assert len(results) == 6
+        assert len(status) == 6
         for status_value in status.values():
-            self.assertEqual(status_value, "optimal")
+            assert status_value == "optimal"
         for reaction, value in results.items():
-            self.assertAlmostEqual(value, expected_results[reaction], 5)
+            assert abs(value - expected_results[reaction]) < 0.00001
 
-    def compare_matrices(self, matrix1, matrix2, places=3):
+    @classmethod
+    def compare_matrices(cls, matrix1, matrix2, places=3):
         nrows = len(matrix1)
         ncols = len(matrix1[0])
-        self.assertEqual(nrows, len(matrix2))
-        self.assertEqual(ncols, len(matrix2[0]))
+        assert nrows == len(matrix2)
+        assert ncols == len(matrix2[0])
         for i in range(nrows):
             for j in range(ncols):
-                self.assertAlmostEqual(matrix1[i][j], matrix2[i][j],
-                                       places=places)
-
-    @skipIf(numpy is None, "double deletions require numpy")
-    def test_double_gene_deletion(self):
-        cobra_model = create_test_model("textbook")
-        genes = ["b0726", "b4025", "b0724", "b0720",
-                 "b2935", "b2935", "b1276", "b1241"]
+                assert abs(matrix1[i][j] - matrix2[i][j]) < 10 ** -places
+
+    @pytest.mark.skipif(numpy is None, reason="double deletions require numpy")
+    def test_double_gene_deletion_benchmark(self, large_model, benchmark):
+        genes = ["b0726", "b4025", "b0724", "b0720", "b2935", "b2935", "b1276",
+                 "b1241"]
+        benchmark(double_gene_deletion, large_model, gene_list1=genes)
+
+    @pytest.mark.skipif(numpy is None, reason="double deletions require numpy")
+    def test_double_gene_deletion(self, model):
+        genes = ["b0726", "b4025", "b0724", "b0720", "b2935", "b2935", "b1276",
+                 "b1241"]
         growth_list = [
             [0.858, 0.857, 0.814, 0.000, 0.858, 0.858, 0.858, 0.858],
             [0.857, 0.863, 0.739, 0.000, 0.863, 0.863, 0.863, 0.863],
@@ -184,73 +188,68 @@ class TestCobraFluxAnalysis(TestCase):
             [0.858, 0.863, 0.814, 0.000, 0.874, 0.874, 0.874, 0.874],
             [0.858, 0.863, 0.814, 0.000, 0.874, 0.874, 0.874, 0.874]]
         opts = {"number_of_processes": 1} if name == "nt" else {}
-        solution = double_gene_deletion(cobra_model, gene_list1=genes, **opts)
-        self.assertEqual(solution["x"], genes)
-        self.assertEqual(solution["y"], genes)
+        solution = double_gene_deletion(model, gene_list1=genes, **opts)
+        assert solution["x"] == genes
+        assert solution["y"] == genes
         self.compare_matrices(growth_list, solution["data"])
         # test when lists differ slightly
-        solution = double_gene_deletion(cobra_model, gene_list1=genes[:-1],
+        solution = double_gene_deletion(model, gene_list1=genes[:-1],
                                         gene_list2=genes,
                                         number_of_processes=1)
-        self.assertEqual(solution["x"], genes[:-1])
-        self.assertEqual(solution["y"], genes)
+        assert solution["x"] == genes[:-1]
+        assert solution["y"] == genes
         self.compare_matrices(growth_list[:-1], solution["data"])
 
-    @skipIf(numpy is None, "double deletions require numpy")
-    def test_double_reaction_deletion(self):
-        cobra_model = create_test_model("textbook")
+    @pytest.mark.skipif(numpy is None, reason="double deletions require numpy")
+    def test_double_reaction_deletion(self, model):
         reactions = ['FBA', 'ATPS4r', 'ENO', 'FRUpts2']
         growth_list = [[0.704, 0.135, 0.000, 0.704],
                        [0.135, 0.374, 0.000, 0.374],
                        [0.000, 0.000, 0.000, 0.000],
                        [0.704, 0.374, 0.000, 0.874]]
 
-        solution = double_reaction_deletion(cobra_model,
+        solution = double_reaction_deletion(model,
                                             reaction_list1=reactions,
                                             number_of_processes=1)
-        self.assertEqual(solution["x"], reactions)
-        self.assertEqual(solution["y"], reactions)
+        assert solution["x"] == reactions
+        assert solution["y"] == reactions
         self.compare_matrices(growth_list, solution["data"])
 
-    def test_flux_variability(self):
-        with open(join(data_directory, "textbook_fva.json"), "r") as infile:
-            fva_results = load(infile)
-
-        infeasible_model = create_test_model("textbook")
+    @pytest.mark.parametrize("solver", list(solver_dict))
+    def test_flux_variability_benchmark(self, large_model, benchmark, solver):
+        benchmark(flux_variability_analysis, large_model, solver=solver,
+                  reaction_list=large_model.reactions[1::3])
+
+    @pytest.mark.parametrize("solver", list(solver_dict))
+    def test_flux_variability(self, model, fva_results, solver):
+        if solver == "esolver":
+            pytest.skip("esolver too slow...")
+        fva_out = flux_variability_analysis(
+            model, solver=solver, reaction_list=model.reactions[1::3])
+        for name, result in iteritems(fva_out):
+            for k, v in iteritems(result):
+                assert abs(fva_results[name][k] - v) < 0.00001
+
+    def test_fva_infeasible(self, model):
+        infeasible_model = model.copy()
         infeasible_model.reactions.get_by_id("EX_glc__D_e").lower_bound = 0
-        for solver in solver_dict:
-            # esolver is really slow
-            if solver == "esolver":
-                continue
-            cobra_model = create_test_model("textbook")
-            fva_out = flux_variability_analysis(
-                cobra_model, solver=solver,
-                reaction_list=cobra_model.reactions[1::3])
-            for name, result in iteritems(fva_out):
-                for k, v in iteritems(result):
-                    self.assertAlmostEqual(fva_results[name][k], v, places=5)
-
-            # ensure that an infeasible model does not run FVA
-            self.assertRaises(ValueError, flux_variability_analysis,
-                              infeasible_model, solver=solver)
-
-    def test_find_blocked_reactions(self):
-        m = create_test_model("textbook")
-        result = find_blocked_reactions(m, m.reactions[40:46])
-        self.assertEqual(result, ['FRUpts2'])
-
-        result = find_blocked_reactions(m, m.reactions[42:48])
-        self.assertEqual(set(result), {'FUMt2_2', 'FRUpts2'})
-
-        result = find_blocked_reactions(m, m.reactions[30:50],
+        # ensure that an infeasible model does not run FVA
+        with pytest.raises(ValueError):
+            flux_variability_analysis(infeasible_model)
+
+    def test_find_blocked_reactions(self, model):
+        result = find_blocked_reactions(model, model.reactions[40:46])
+        assert result == ['FRUpts2']
+
+        result = find_blocked_reactions(model, model.reactions[42:48])
+        assert set(result) == {'FUMt2_2', 'FRUpts2'}
+
+        result = find_blocked_reactions(model, model.reactions[30:50],
                                         open_exchanges=True)
-        self.assertEqual(result, [])
+        assert result == []
 
-    def test_loopless(self):
-        try:
-            solver = get_solver_name(mip=True)
-        except:
-            self.skipTest("no MILP solver found")
+    @classmethod
+    def construct_ll_test_model(cls):
         test_model = Model()
         test_model.add_metabolites(Metabolite("A"))
         test_model.add_metabolites(Metabolite("B"))
@@ -270,17 +269,29 @@ class TestCobraFluxAnalysis(TestCase):
                             test_model.metabolites.A: 1})
         DM_C.objective_coefficient = 1
         test_model.add_reactions([EX_A, DM_C, v1, v2, v3])
+        return test_model
+
+    def test_loopless_benchmark(self, benchmark):
+        test_model = self.construct_ll_test_model()
+        benchmark(lambda: construct_loopless_model(test_model).optimize())
+
+    def test_loopless(self):
+        try:
+            get_solver_name(mip=True)
+        except SolverNotFound:
+            pytest.skip("no MILP solver found")
+        test_model = self.construct_ll_test_model()
         feasible_sol = construct_loopless_model(test_model).optimize()
-        v3.lower_bound = 1
+        test_model.reactions.get_by_id('v3').lower_bound = 1
         infeasible_sol = construct_loopless_model(test_model).optimize()
-        self.assertEqual(feasible_sol.status, "optimal")
-        self.assertEqual(infeasible_sol.status, "infeasible")
+        assert feasible_sol.status == "optimal"
+        assert infeasible_sol.status == "infeasible"
 
     def test_gapfilling(self):
         try:
-            solver = get_solver_name(mip=True)
-        except:
-            self.skipTest("no MILP solver found")
+            get_solver_name(mip=True)
+        except SolverNotFound:
+            pytest.skip("no MILP solver found")
         m = Model()
         m.add_metabolites(map(Metabolite, ["a", "b", "c"]))
         r = Reaction("EX_A")
@@ -304,56 +315,52 @@ class TestCobraFluxAnalysis(TestCase):
 
         # GrowMatch
         result = gapfilling.growMatch(m, U)[0]
-        self.assertEqual(len(result), 1)
-        self.assertEqual(result[0].id, "a2b")
+        assert len(result) == 1
+        assert result[0].id == "a2b"
         # SMILEY
         result = gapfilling.SMILEY(m, "b", U)[0]
-        self.assertEqual(len(result), 1)
-        self.assertEqual(result[0].id, "a2b")
+        assert len(result) == 1
+        assert result[0].id == "a2b"
 
         # 2 rounds of GrowMatch with exchange reactions
         result = gapfilling.growMatch(m, None, ex_rxns=True, iterations=2)
-        self.assertEqual(len(result), 2)
-        self.assertEqual(len(result[0]), 1)
-        self.assertEqual(len(result[1]), 1)
-        self.assertEqual({i[0].id for i in result},
-                         {"SMILEY_EX_b", "SMILEY_EX_c"})
-
-    @skipIf(numpy is None, "phase plane requires numpy")
-    def test_phenotype_phase_plane(self):
-        model = create_test_model("textbook")
+        assert len(result) == 2
+        assert len(result[0]) == 1
+        assert len(result[1]) == 1
+        assert {i[0].id for i in result} == {"SMILEY_EX_b", "SMILEY_EX_c"}
+
+    @pytest.mark.skipif(numpy is None, reason="phase plane require numpy")
+    def test_phenotype_phase_plane_benchmark(self, model, benchmark):
+        benchmark(calculate_phenotype_phase_plane,
+                  model, "EX_glc__D_e", "EX_o2_e",
+                  reaction1_npoints=20, reaction2_npoints=20)
+
+    @pytest.mark.skipif(numpy is None, reason="phase plane require numpy")
+    def test_phenotype_phase_plane(self, model):
         data = calculate_phenotype_phase_plane(
             model, "EX_glc__D_e", "EX_o2_e",
             reaction1_npoints=20, reaction2_npoints=20)
-        self.assertEqual(data.growth_rates.shape, (20, 20))
-        self.assertAlmostEqual(data.growth_rates.max(), 1.20898, places=4)
-        self.assertAlmostEqual(abs(data.growth_rates[0, :]).max(), 0, places=4)
+        assert data.growth_rates.shape == (20, 20)
+        assert abs(data.growth_rates.max() - 1.20898) < 0.0001
+        assert abs(data.growth_rates[0, :].max()) < 0.0001
         if matplotlib is None:
-            self.skipTest("can't test plots without matplotlib")
+            pytest.skip("can't test plots without matplotlib")
         data.plot()
 
     def check_entries(self, out, desired_entries):
         """ensure each entry in desired_entries appears in output"""
         output = out.getvalue().strip()
         output_set = set((re.sub('\s', '', l) for l in output.splitlines()))
-
         for item in desired_entries:
-            self.assertIn(re.sub('\s', '', item), output_set)
-
-    @skipIf((pandas is None) or (tabulate is None),
-            "summary methods require pandas and tabulate")
-    def test_summary_methods(self):
+            assert re.sub('\s', '', item) in output_set
 
+    @pytest.mark.skipif((pandas is None) or (tabulate is None),
+                        reason="summary methods require pandas and tabulate")
+    def test_summary_methods(self, model, solved_model):
         # Test model summary methods
-        model = create_test_model("textbook")
-        with self.assertRaises(Exception):
+        with pytest.raises(Exception):
             model.summary()
 
-        # Load model solution
-        with open(join(data_directory, "textbook_solution.pickle"),
-                  "rb") as infile:
-            model.solution = pickle.load(infile)
-
         desired_entries = [
             'idFluxRangeidFluxRangeBiomass_Ecol...0.874',
             'o2_e       21.8   [19.9, 23.7]'
@@ -375,7 +382,7 @@ class TestCobraFluxAnalysis(TestCase):
         ]
         for solver in solver_dict:
             with captured_output() as (out, err):
-                model.summary(fva=0.95, solver=solver)
+                solved_model.summary(fva=0.95, solver=solver)
             self.check_entries(out, desired_entries)
 
         # test non-fva version (these should be fixed for textbook model
@@ -393,11 +400,11 @@ class TestCobraFluxAnalysis(TestCase):
         # there are multiple entries per line.
         for solver in solver_dict:
             with captured_output() as (out, err):
-                model.summary()
+                solved_model.summary()
 
             s = out.getvalue()
             for i in desired_entries:
-                self.assertIn(i, s)
+                assert i in s
 
         # Test metabolite summary methods
         desired_entries = [
@@ -414,7 +421,7 @@ class TestCobraFluxAnalysis(TestCase):
 
         for solver in solver_dict:
             with captured_output() as (out, err):
-                model.metabolites.q8_c.summary()
+                solved_model.metabolites.q8_c.summary()
             self.check_entries(out, desired_entries)
 
         desired_entries = [
@@ -433,17 +440,5 @@ class TestCobraFluxAnalysis(TestCase):
 
         for solver in solver_dict:
             with captured_output() as (out, err):
-                model.metabolites.fdp_c.summary(fva=0.99, solver=solver)
+                solved_model.metabolites.fdp_c.summary(fva=0.99, solver=solver)
             self.check_entries(out, desired_entries)
-
-
-# make a test suite to run all of the tests
-loader = TestLoader()
-suite = loader.loadTestsFromModule(sys.modules[__name__])
-
-
-def test_all():
-    TextTestRunner(verbosity=2).run(suite)
-
-if __name__ == "__main__":
-    test_all()
diff --git a/cobra/test/test_io.py b/cobra/test/test_io.py
new file mode 100644
index 0000000..4649a52
--- /dev/null
+++ b/cobra/test/test_io.py
@@ -0,0 +1,213 @@
+from warnings import warn
+from tempfile import gettempdir
+from os import unlink
+from os.path import join, split
+from six import iteritems
+from functools import partial
+from pickle import load, dump
+import pytest
+from collections import namedtuple
+from cobra import io
+from .conftest import data_directory
+
+
+def write_legacy_sbml_placeholder():
+    pass
+
+
+try:
+    import scipy
+except ImportError:
+    scipy = None
+try:
+    import libsbml
+
+    write_legacy_sbml = io.write_legacy_sbml
+except ImportError:
+    libsbml = None
+    write_legacy_sbml = write_legacy_sbml_placeholder
+try:
+    import jsonschema
+except ImportError:
+    jsonschema = None
+try:
+    import cPickle
+
+    cload = cPickle.load
+    cdump = cPickle.dump
+except ImportError:
+    cPickle = None
+    cload = None
+    cdump = None
+
+
+def validate_json(filename):
+    with open(filename, "r") as infile:
+        loaded = io.json.json.load(infile)
+    if jsonschema is None:
+        warn("jsonschema not installed")
+    else:
+        jsonschema.validate(loaded, io.json.json_schema)
+
+
+def read_pickle(filename, load_function=load):
+    with open(filename, "rb") as infile:
+        return load_function(infile)
+
+
+def write_pickle(model, filename, dump_function=dump):
+    with open(filename, "wb") as outfile:
+        dump_function(model, outfile)
+
+
+IOTrial = namedtuple('IOTrial',
+                     ['name', 'reference_file', 'test_file', 'read_function',
+                      'write_function', 'validation_function'])
+trials = [IOTrial('fbc2', 'mini.pickle', 'mini_fbc2.xml',
+                  io.read_sbml_model, io.write_sbml_model,
+                  io.sbml3.validate_sbml_model),
+          IOTrial('fbc2Gz', 'mini.pickle', 'mini_fbc2.xml.gz',
+                  io.read_sbml_model, io.write_sbml_model, None),
+          IOTrial('fbc2Bz2', 'mini.pickle', 'mini_fbc2.xml.bz2',
+                  io.read_sbml_model, io.write_sbml_model, None),
+          pytest.mark.skipif("not libsbml")(
+              IOTrial('fbc1', 'mini.pickle', 'mini_fbc1.xml',
+                      io.read_sbml_model,
+                      partial(write_legacy_sbml, use_fbc_package=True), None)),
+          pytest.mark.skipif("not libsbml")(
+              IOTrial('cobra', 'mini.pickle', 'mini_cobra.xml',
+                      io.read_sbml_model,
+                      partial(write_legacy_sbml, use_fbc_package=False),
+                      None)),
+          pytest.mark.skipif("not scipy")(
+              IOTrial('mat', 'mini.pickle', 'mini.mat',
+                      io.load_matlab_model, io.save_matlab_model, None)),
+          IOTrial('json', 'mini.pickle', 'mini.json',
+                  io.load_json_model, io.save_json_model, validate_json),
+          IOTrial('pickle', 'mini.pickle', 'mini.pickle',
+                  read_pickle, write_pickle, None),
+          pytest.mark.skipif("not cPickle")(
+              IOTrial('pickle', 'mini.pickle', 'mini.pickle',
+                      partial(read_pickle, load_function=cload),
+                      partial(write_pickle, dump_function=cdump), None))
+          ]
+trial_names = [node.name for node in trials]
+
+
+ at pytest.fixture(scope="module", params=trials, ids=trial_names)
+def io_trial(request, data_directory):
+    with open(join(data_directory, request.param.reference_file),
+              "rb") as infile:
+        reference_model = load(infile)
+    test_model = request.param.read_function(join(data_directory,
+                                                  request.param.test_file))
+    test_output_filename = join(gettempdir(),
+                                split(request.param.test_file)[-1])
+    request.param.write_function(test_model, test_output_filename)
+    reread_model = request.param.read_function(test_output_filename)
+    unlink(test_output_filename)
+    return request.param.name, reference_model, test_model, reread_model
+
+
+class TestCobraIO:
+    @classmethod
+    def compare_models(cls, name, model1, model2):
+        assert len(model1.reactions) == len(model2.reactions)
+        assert len(model1.metabolites) == len(model2.metabolites)
+        for attr in ("id", "name", "lower_bound", "upper_bound",
+                     "objective_coefficient", "gene_reaction_rule"):
+            assert getattr(model1.reactions[0], attr) == getattr(
+                model2.reactions[0], attr)
+            assert getattr(model1.reactions[10], attr) == getattr(
+                model2.reactions[10], attr)
+            assert getattr(model1.reactions[-1], attr) == getattr(
+                model2.reactions[-1], attr)
+        for attr in ("id", "name", "compartment", "formula", "charge"):
+            assert getattr(model1.metabolites[0], attr) == getattr(
+                model2.metabolites[0], attr)
+            assert getattr(model1.metabolites[10], attr) == getattr(
+                model2.metabolites[10], attr)
+            assert getattr(model1.metabolites[-1], attr) == getattr(
+                model2.metabolites[-1], attr)
+        assert len(model1.reactions[0].metabolites) == len(
+            model2.reactions[0].metabolites)
+        assert len(model1.reactions[14].metabolites) == len(
+            model2.reactions[14].metabolites)
+        assert len(model1.reactions[-1].metabolites) == len(
+            model2.reactions[-1].metabolites)
+        assert len(model1.genes) == len(model2.genes)
+        # ensure they have the same solution max
+        model1.optimize()
+        model2.optimize()
+        assert abs(model1.solution.f - model2.solution.f) < 0.001
+        # ensure the references are correct
+        if name == 'mat':
+            pytest.xfail('curiously, reference check now fail with matlab')
+        assert model2.metabolites[0]._model is model2
+        assert model2.reactions[0]._model is model2
+        assert model2.genes[0]._model is model2
+
+    @classmethod
+    def extra_comparisons(cls, name, model1, model2):
+        if name in ['fbc1', 'mat', 'cobra']:
+            pytest.skip("not supported")
+        assert model1.compartments == model2.compartments
+        assert model1.metabolites[4].annotation == model2.metabolites[
+            4].annotation
+        assert model1.reactions[4].annotation == model2.reactions[4].annotation
+        assert model1.genes[4].annotation == model2.genes[4].annotation
+        for attr in ("id", "name"):
+            assert getattr(model1.genes[0], attr) == getattr(model2.genes[0],
+                                                             attr)
+            assert getattr(model1.genes[10], attr) == getattr(model2.genes[10],
+                                                              attr)
+            assert getattr(model1.genes[-1], attr) == getattr(model2.genes[-1],
+                                                              attr)
+
+    def test_read(self, io_trial):
+        name, reference_model, test_model, _ = io_trial
+        if name in ['fbc1']:
+            pytest.xfail('not supported')
+        self.compare_models(name, reference_model, test_model)
+        self.extra_comparisons(name, reference_model, test_model)
+
+    def test_write(self, io_trial):
+        name, _, test_model, reread_model = io_trial
+        if name in ['fbc1']:
+            pytest.xfail('not supported')
+        self.compare_models(name, test_model, reread_model)
+        self.extra_comparisons(name, test_model, reread_model)
+
+
+def test_benchmark_read(data_directory, benchmark):
+    benchmark(io.sbml3.read_sbml_model, join(data_directory, 'mini_fbc2.xml'))
+
+
+def test_benchmark_write(model, benchmark):
+    benchmark(io.sbml3.write_sbml_model, model, join(gettempdir(), "-bench"))
+
+
+ at pytest.mark.parametrize("trial", trials)
+def test_validate(trial, data_directory):
+    if trial.validation_function is None:
+        pytest.skip('not implemented')
+    test_file = join(data_directory, trial.test_file)
+    trial.validation_function(test_file)
+
+
+ at pytest.mark.parametrize("trial", trials)
+def test_read_nonexistent(trial):
+    pytest.raises(IOError, trial.read_function, "fake_file")
+
+
+def test_sbml_error(data_directory):
+    filename = join(data_directory, "invalid0.xml")
+    with pytest.raises(io.sbml3.CobraSBMLError):
+        io.read_sbml_model(filename)
+
+
+def test_bad_validation(data_directory):
+    for i in range(3):
+        filename = join(data_directory, "invalid%d.xml" % i)
+        m, errors = io.sbml3.validate_sbml_model(filename)
+        assert any(len(v) >= 1 for v in iteritems(errors))
diff --git a/cobra/test/manipulation.py b/cobra/test/test_manipulation.py
similarity index 52%
rename from cobra/test/manipulation.py
rename to cobra/test/test_manipulation.py
index 30f5c08..6bafd08 100644
--- a/cobra/test/manipulation.py
+++ b/cobra/test/test_manipulation.py
@@ -1,115 +1,95 @@
-from unittest import TestCase, TestLoader, TextTestRunner
+from cobra.core import Metabolite, Model, Reaction
+from cobra.manipulation import *
+import pytest
+from .conftest import model, salmonella
 
-import sys
 
-if __name__ == "__main__":
-    sys.path.insert(0, "../..")
-    from cobra.test import create_test_model, data_directory
-    from cobra.core import Metabolite, Model, Reaction
-    from cobra.manipulation import *
-    sys.path.pop(0)
-else:
-    from . import create_test_model, data_directory
-    from ..core import Metabolite, Model, Reaction
-    from ..manipulation import *
-
-
-class TestManipulation(TestCase):
+class TestManipulation:
     """Test functions in cobra.manipulation"""
 
-    def test_canonical_form(self):
-        model = create_test_model("textbook")
+    def test_canonical_form(self, model):
         # add G constraint to test
         g_constr = Metabolite("SUCCt2_2__test_G_constraint")
         g_constr._constraint_sense = "G"
         g_constr._bound = 5.0
         model.reactions.get_by_id("SUCCt2_2").add_metabolites({g_constr: 1})
-        self.assertAlmostEqual(model.optimize("maximize").f, 0.855, places=3)
+        assert abs(model.optimize("maximize").f - 0.855) < 0.001
         # convert to canonical form
         model = canonical_form(model)
-        self.assertAlmostEqual(model.optimize("maximize").f, 0.855, places=3)
+        assert abs(model.optimize("maximize").f - 0.855) < 10 ** -3
 
-    def test_canonical_form_minimize(self):
-        model = create_test_model("textbook")
+    def test_canonical_form_minimize(self, model):
         # make a minimization problem
         model.reactions.get_by_id("Biomass_Ecoli_core").lower_bound = 0.5
         for reaction in model.reactions:
             reaction.objective_coefficient = reaction.id == "GAPD"
-        self.assertAlmostEqual(model.optimize("minimize").f, 6.27, places=3)
+        assert abs(model.optimize("minimize").f - 6.27) < 10 ** -3
         # convert to canonical form. Convert minimize to maximize
         model = canonical_form(model, objective_sense="minimize")
-        self.assertAlmostEqual(model.optimize("maximize").f, -6.27, places=3)
+        assert abs(model.optimize("maximize").f + 6.27) < 10 ** -3
         # lower bounds should now be <= constraints
-        self.assertEqual(
-            model.reactions.get_by_id("Biomass_Ecoli_core").lower_bound, 0.0)
+        assert model.reactions.get_by_id(
+            "Biomass_Ecoli_core").lower_bound == 0.0
 
-    def test_modify_reversible(self):
-        model1 = create_test_model("textbook")
+    def test_modify_reversible(self, model):
+        model1 = model.copy()
         model1.optimize()
-        model2 = create_test_model("textbook")
+        model2 = model.copy()
         convert_to_irreversible(model2)
         model2.optimize()
-        self.assertAlmostEqual(model1.solution.f, model2.solution.f, places=3)
+        assert abs(model1.solution.f - model2.solution.f) < 10 ** -3
         revert_to_reversible(model2)
         model2.optimize()
-        self.assertAlmostEqual(model1.solution.f, model2.solution.f, places=3)
-
+        assert abs(model1.solution.f - model2.solution.f) < 10 ** -3
         # Ensure revert_to_reversible is robust to solutions generated both
         # before and after reversibility conversion, or not solved at all.
-        model3 = create_test_model("textbook")
+        model3 = model.copy()
         model3.optimize()
         convert_to_irreversible(model3)
         revert_to_reversible(model3)
-        self.assertAlmostEqual(model1.solution.f, model3.solution.f, places=3)
-
+        assert abs(model1.solution.f - model3.solution.f) < 10 ** -3
         # test reaction where both bounds are negative
-        model4 = create_test_model("textbook")
+        model4 = model.copy()
         glc = model4.reactions.get_by_id("EX_glc__D_e")
         glc.upper_bound = -1
         convert_to_irreversible(model4)
         model4.optimize()
-        self.assertAlmostEqual(model1.solution.f, model4.solution.f, places=3)
+        assert abs(model1.solution.f - model4.solution.f) < 10 ** -3
         glc_rev = model4.reactions.get_by_id(glc.notes["reflection"])
-        self.assertEqual(glc_rev.lower_bound, 1)
-        self.assertEqual(glc.upper_bound, 0)
+        assert glc_rev.lower_bound == 1
+        assert glc.upper_bound == 0
         revert_to_reversible(model4)
-        self.assertEqual(glc.upper_bound, -1)
+        assert glc.upper_bound == -1
 
-    def test_escape_ids(self):
-        model = create_test_model('textbook')
+    def test_escape_ids(self, model):
         model.reactions.PGI.gene_reaction_rule = "a.b or c"
-        self.assertIn("a.b", model.genes)
+        assert "a.b" in model.genes
         escape_ID(model)
-        self.assertNotIn("a.b", model.genes)
+        assert "a.b" not in model.genes
 
-    def test_rename_gene(self):
-        model = create_test_model('textbook')
+    def test_rename_gene(self, model):
         original_name = model.genes.b1241.name
         rename_dict = {"b1241": "foo", "hello": "world",
                        "b2465": "b3919", "bar": "2935"}
         modify.rename_genes(model, rename_dict)
         for i in rename_dict:
-            self.assertNotIn(i, model.genes)
-        self.assertIn("foo", model.genes)
+            assert i not in model.genes
+        assert "foo" in model.genes
         # make sure the object name was preserved
-        self.assertEqual(model.genes.foo.name, original_name)
+        assert model.genes.foo.name == original_name
         # make sure the reactions are correct
-        self.assertEqual(len(model.genes.foo.reactions), 2)
-        self.assertEqual(model.reactions.ACALD.gene_reaction_rule,
-                         "b0351 or foo")
-        self.assertEqual(model.reactions.TPI.gene_reaction_rule, "b3919")
-        self.assertEqual(model.reactions.TPI.genes, {model.genes.b3919})
-        self.assertEqual(model.reactions.TKT1.gene_reaction_rule,
-                         "b2935 or b3919")
-        self.assertEqual(model.reactions.TKT1.genes,
-                         {model.genes.b2935, model.genes.b3919})
-        self.assertEqual(model.genes.b3919.reactions,
-                         {model.reactions.get_by_id(i)
-                          for i in ("TKT1", "TKT2", "TPI")})
-
-    def test_gene_knockout_computation(self):
-        cobra_model = create_test_model()
-
+        assert len(model.genes.foo.reactions) == 2
+        assert model.reactions.ACALD.gene_reaction_rule == "b0351 or foo"
+        assert model.reactions.TPI.gene_reaction_rule == "b3919"
+        assert model.reactions.TPI.genes == {model.genes.b3919}
+        assert model.reactions.TKT1.gene_reaction_rule == "b2935 or b3919"
+        assert model.reactions.TKT1.genes == {
+            model.genes.b2935, model.genes.b3919}
+        assert model.genes.b3919.reactions == {model.reactions.get_by_id(i)
+                                               for i in
+                                               ("TKT1", "TKT2", "TPI")}
+
+    def test_gene_knockout_computation(self, salmonella):
         def find_gene_knockout_reactions_fast(cobra_model, gene_list):
             compiled_rules = get_compiled_gene_reaction_rules(
                 cobra_model)
@@ -126,10 +106,10 @@ class TestManipulation(TestCase):
                                   for i in expected_reaction_ids}
             removed1 = set(find_gene_knockout_reactions(m, genes))
             removed2 = set(find_gene_knockout_reactions_fast(m, genes))
-            self.assertEqual(removed1, expected_reactions)
-            self.assertEqual(removed2, expected_reactions)
+            assert removed1 == expected_reactions
+            assert removed2 == expected_reactions
             delete_model_genes(m, gene_ids, cumulative_deletions=False)
-            self.assertEqual(get_removed(m), expected_reaction_ids)
+            assert get_removed(m) == expected_reaction_ids
             undelete_model_genes(m)
 
         gene_list = ['STM1067', 'STM0227']
@@ -137,24 +117,24 @@ class TestManipulation(TestCase):
                                '3HAD180', '3HAD100', '3HAD181', '3HAD120',
                                '3HAD60', '3HAD141', '3HAD161', 'T2DECAI',
                                '3HAD40'}
-        test_computation(cobra_model, gene_list, dependent_reactions)
-        test_computation(cobra_model, ['STM4221'], {'PGI'})
-        test_computation(cobra_model, ['STM1746.S'], {'4PEPTabcpp'})
+        test_computation(salmonella, gene_list, dependent_reactions)
+        test_computation(salmonella, ['STM4221'], {'PGI'})
+        test_computation(salmonella, ['STM1746.S'], {'4PEPTabcpp'})
         # test cumulative behavior
-        delete_model_genes(cobra_model, gene_list[:1])
-        delete_model_genes(cobra_model, gene_list[1:],
+        delete_model_genes(salmonella, gene_list[:1])
+        delete_model_genes(salmonella, gene_list[1:],
                            cumulative_deletions=True)
-        delete_model_genes(cobra_model, ["STM4221"],
+        delete_model_genes(salmonella, ["STM4221"],
                            cumulative_deletions=True)
         dependent_reactions.add('PGI')
-        self.assertEqual(get_removed(cobra_model), dependent_reactions)
+        assert get_removed(salmonella) == dependent_reactions
         # non-cumulative following cumulative
-        delete_model_genes(cobra_model, ["STM4221"],
+        delete_model_genes(salmonella, ["STM4221"],
                            cumulative_deletions=False)
-        self.assertEqual(get_removed(cobra_model), {'PGI'})
+        assert get_removed(salmonella) == {'PGI'}
         # make sure on reset that the bounds are correct
-        reset_bound = cobra_model.reactions.get_by_id("T2DECAI").upper_bound
-        self.assertEqual(reset_bound, 1000.)
+        reset_bound = salmonella.reactions.get_by_id("T2DECAI").upper_bound
+        assert reset_bound == 1000.
         # test computation when gene name is a subset of another
         test_model = Model()
         test_reaction_1 = Reaction("test1")
@@ -183,7 +163,7 @@ class TestManipulation(TestCase):
     def test_remove_genes(self):
         m = Model("test")
         m.add_reactions([Reaction("r" + str(i + 1)) for i in range(8)])
-        self.assertEqual(len(m.reactions), 8)
+        assert len(m.reactions) == 8
         rxns = m.reactions
         rxns.r1.gene_reaction_rule = "(a and b) or (c and a)"
         rxns.r2.gene_reaction_rule = "(a and b and d and e)"
@@ -193,34 +173,33 @@ class TestManipulation(TestCase):
         rxns.r6.gene_reaction_rule = "y"
         rxns.r7.gene_reaction_rule = "x or     z"
         rxns.r8.gene_reaction_rule = ""
-        self.assertIn("a", m.genes)
-        self.assertIn("x", m.genes)
+        assert "a" in m.genes
+        assert "x" in m.genes
         remove_genes(m, ["a"], remove_reactions=False)
-        self.assertNotIn("a", m.genes)
-        self.assertIn("x", m.genes)
-        self.assertEqual(rxns.r1.gene_reaction_rule, "")
-        self.assertEqual(rxns.r2.gene_reaction_rule, "")
-        self.assertEqual(rxns.r3.gene_reaction_rule, "b and c")
-        self.assertEqual(rxns.r4.gene_reaction_rule, "(f and b) or (b and c)")
-        self.assertEqual(rxns.r5.gene_reaction_rule, "x")
-        self.assertEqual(rxns.r6.gene_reaction_rule, "y")
-        self.assertEqual(rxns.r7.genes, {m.genes.x, m.genes.z})
-        self.assertEqual(rxns.r8.gene_reaction_rule, "")
+        assert "a" not in m.genes
+        assert "x" in m.genes
+        assert rxns.r1.gene_reaction_rule == ""
+        assert rxns.r2.gene_reaction_rule == ""
+        assert rxns.r3.gene_reaction_rule == "b and c"
+        assert rxns.r4.gene_reaction_rule == "(f and b) or (b and c)"
+        assert rxns.r5.gene_reaction_rule == "x"
+        assert rxns.r6.gene_reaction_rule == "y"
+        assert rxns.r7.genes == {m.genes.x, m.genes.z}
+        assert rxns.r8.gene_reaction_rule == ""
         remove_genes(m, ["x"], remove_reactions=True)
-        self.assertEqual(len(m.reactions), 7)
-        self.assertNotIn("r5", m.reactions)
-        self.assertNotIn("x", m.genes)
-        self.assertEqual(rxns.r1.gene_reaction_rule, "")
-        self.assertEqual(rxns.r2.gene_reaction_rule, "")
-        self.assertEqual(rxns.r3.gene_reaction_rule, "b and c")
-        self.assertEqual(rxns.r4.gene_reaction_rule, "(f and b) or (b and c)")
-        self.assertEqual(rxns.r6.gene_reaction_rule, "y")
-        self.assertEqual(rxns.r7.gene_reaction_rule, "z")
-        self.assertEqual(rxns.r7.genes, {m.genes.z})
-        self.assertEqual(rxns.r8.gene_reaction_rule, "")
-
-    def test_SBO_annotation(self):
-        model = create_test_model("textbook")
+        assert len(m.reactions) == 7
+        assert "r5" not in m.reactions
+        assert "x" not in m.genes
+        assert rxns.r1.gene_reaction_rule == ""
+        assert rxns.r2.gene_reaction_rule == ""
+        assert rxns.r3.gene_reaction_rule == "b and c"
+        assert rxns.r4.gene_reaction_rule == "(f and b) or (b and c)"
+        assert rxns.r6.gene_reaction_rule == "y"
+        assert rxns.r7.gene_reaction_rule == "z"
+        assert rxns.r7.genes == {m.genes.z}
+        assert rxns.r8.gene_reaction_rule == ""
+
+    def test_sbo_annotation(self, model):
         rxns = model.reactions
         rxns.EX_o2_e.annotation.clear()
         fake_DM = Reaction("DM_h_c")
@@ -230,46 +209,36 @@ class TestManipulation(TestCase):
         # an existing SBO annotation
         rxns.get_by_id("EX_h_e").annotation["SBO"] = "SBO:0000628"
         add_SBO(model)
-        self.assertEqual(rxns.EX_o2_e.annotation["SBO"], "SBO:0000627")
-        self.assertEqual(rxns.DM_h_c.annotation["SBO"], "SBO:0000628")
-        self.assertEqual(rxns.EX_h_e.annotation["SBO"], "SBO:0000628")
+        assert rxns.EX_o2_e.annotation["SBO"] == "SBO:0000627"
+        assert rxns.DM_h_c.annotation["SBO"] == "SBO:0000628"
+        assert rxns.EX_h_e.annotation["SBO"] == "SBO:0000628"
 
-    def test_validate_reaction_bounds(self):
-        model = create_test_model("textbook")
+    def test_validate_reaction_bounds(self, model):
         model.reactions[0].lower_bound = float("-inf")
         model.reactions[1].lower_bound = float("nan")
         model.reactions[0].upper_bound = float("inf")
         model.reactions[1].upper_bound = float("nan")
         errors = check_reaction_bounds(model)
-        self.assertEqual(len(errors), 4)
+        assert len(errors) == 4
 
-    def test_validate_formula_compartment(self):
-        model = create_test_model("textbook")
+    def test_validate_formula_compartment(self, model):
         model.metabolites[1].compartment = "fake"
         model.metabolites[1].formula = "(a*.bcde)"
         errors = check_metabolite_compartment_formula(model)
-        self.assertEqual(len(errors), 2)
+        assert len(errors) == 2
 
-    def test_validate_mass_balance(self):
-        model = create_test_model("textbook")
-        self.assertEqual(len(check_mass_balance(model)), 0)
+    def test_validate_mass_balance(self, model):
+        assert len(check_mass_balance(model)) == 0
         # if we remove the SBO term which marks the reaction as
         # mass balanced, then the reaction should be detected as
         # no longer mass balanced
         EX_rxn = model.reactions.query("EX")[0]
         EX_rxn.annotation.pop("SBO")
         balance = check_mass_balance(model)
-        self.assertEqual(len(balance), 1)
-        self.assertIn(EX_rxn, balance)
-
-
-# make a test suite to run all of the tests
-loader = TestLoader()
-suite = loader.loadTestsFromModule(sys.modules[__name__])
-
-
-def test_all():
-    TextTestRunner(verbosity=2).run(suite)
-
-if __name__ == "__main__":
-    test_all()
+        assert len(balance) == 1
+        assert EX_rxn in balance
+        m1 = Metabolite('m1', formula='()')
+        r1 = Reaction('r1')
+        r1.add_metabolites({m1: 1})
+        with pytest.raises(ValueError):
+            r1.check_mass_balance()
diff --git a/cobra/test/test_model.py b/cobra/test/test_model.py
new file mode 100644
index 0000000..56d8904
--- /dev/null
+++ b/cobra/test/test_model.py
@@ -0,0 +1,597 @@
+from copy import deepcopy
+import warnings
+import pytest
+from cobra.core import Model, Metabolite, Reaction
+from cobra.solvers import solver_dict
+from .conftest import model, array_model
+
+try:
+    import scipy
+except ImportError:
+    scipy = None
+
+
+class TestReactions:
+    def test_gpr(self):
+        model = Model()
+        reaction = Reaction("test")
+        # set a gpr to  reaction not in a model
+        reaction.gene_reaction_rule = "(g1 or g2) and g3"
+        assert reaction.gene_reaction_rule == "(g1 or g2) and g3"
+        assert len(reaction.genes) == 3
+        # adding reaction with a GPR propagates to the model
+        model.add_reaction(reaction)
+        assert len(model.genes) == 3
+        # ensure the gene objects are the same in the model and reaction
+        reaction_gene = list(reaction.genes)[0]
+        model_gene = model.genes.get_by_id(reaction_gene.id)
+        assert reaction_gene is model_gene
+        # test ability to handle uppercase AND/OR
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore")
+            reaction.gene_reaction_rule = "(b1 AND b2) OR (b3 and b4)"
+        assert reaction.gene_reaction_rule == "(b1 and b2) or (b3 and b4)"
+        assert len(reaction.genes) == 4
+        # ensure regular expressions correctly extract genes from malformed
+        # GPR string
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore")
+            reaction.gene_reaction_rule = "(a1 or a2"
+            assert len(reaction.genes) == 2
+            reaction.gene_reaction_rule = "(forT or "
+            assert len(reaction.genes) == 1
+
+    def test_gpr_modification(self, model):
+        reaction = model.reactions.get_by_id("PGI")
+        old_gene = list(reaction.genes)[0]
+        new_gene = model.genes.get_by_id("s0001")
+        # add an existing 'gene' to the gpr
+        reaction.gene_reaction_rule = 's0001'
+        assert new_gene in reaction.genes
+        assert reaction in new_gene.reactions
+        # removed old gene correctly
+        assert old_gene not in reaction.genes
+        assert reaction not in old_gene.reactions
+        # add a new 'gene' to the gpr
+        reaction.gene_reaction_rule = 'fake_gene'
+        assert model.genes.has_id("fake_gene")
+        fake_gene = model.genes.get_by_id("fake_gene")
+        assert fake_gene in reaction.genes
+        assert reaction in fake_gene.reactions
+        fake_gene.name = "foo_gene"
+        assert reaction.gene_name_reaction_rule == fake_gene.name
+
+    @pytest.mark.parametrize("solver", list(solver_dict))
+    def test_add_metabolite_benchmark(self, model, benchmark, solver):
+        reaction = model.reactions.get_by_id("PGI")
+        many_metabolites = dict((m, 1) for m in model.metabolites[0:50])
+
+        def add_remove_metabolite():
+            reaction.add_metabolites(many_metabolites)
+            if not getattr(model, 'solver', None):
+                solver_dict[solver].create_problem(model)
+            for m, c in many_metabolites.items():
+                try:
+                    reaction.pop(m.id)
+                except KeyError:
+                    pass
+        benchmark(add_remove_metabolite)
+
+    def test_add_metabolite(self, model):
+        reaction = model.reactions.get_by_id("PGI")
+        reaction.add_metabolites({model.metabolites[0]: 1})
+        assert model.metabolites[0] in reaction._metabolites
+        fake_metabolite = Metabolite("fake")
+        reaction.add_metabolites({fake_metabolite: 1})
+        assert fake_metabolite in reaction._metabolites
+        assert model.metabolites.has_id("fake")
+        assert model.metabolites.get_by_id("fake") is fake_metabolite
+
+        # test adding by string
+        reaction.add_metabolites({"g6p_c": -1})  # already in reaction
+        assert reaction._metabolites[
+                   model.metabolites.get_by_id("g6p_c")] == -2
+        reaction.add_metabolites({"h_c": 1})
+        assert reaction._metabolites[model.metabolites.get_by_id("h_c")] == 1
+        with pytest.raises(KeyError):
+            reaction.add_metabolites({"missing": 1})
+
+        # test adding to a new Reaction
+        reaction = Reaction("test")
+        assert len(reaction._metabolites) == 0
+        reaction.add_metabolites({Metabolite("test_met"): -1})
+        assert len(reaction._metabolites) == 1
+
+    @pytest.mark.parametrize("solver", list(solver_dict))
+    def test_subtract_metabolite_benchmark(self, model, benchmark, solver):
+        benchmark(self.test_subtract_metabolite, model, solver)
+
+    @pytest.mark.parametrize("solver", list(solver_dict))
+    def test_subtract_metabolite(self, model, solver):
+        reaction = model.reactions.get_by_id("PGI")
+        reaction.subtract_metabolites(reaction.metabolites)
+        if not getattr(model, 'solver', None):
+            solver_dict[solver].create_problem(model)
+        assert len(reaction.metabolites) == 0
+
+    def test_mass_balance(self, model):
+        reaction = model.reactions.get_by_id("PGI")
+        # should be balanced now
+        assert len(reaction.check_mass_balance()) == 0
+        # should not be balanced after adding a hydrogen
+        reaction.add_metabolites({model.metabolites.get_by_id("h_c"): 1})
+        imbalance = reaction.check_mass_balance()
+        assert "H" in imbalance
+        assert imbalance["H"] == 1
+
+    def test_build_from_string(self, model):
+        m = len(model.metabolites)
+        pgi = model.reactions.get_by_id("PGI")
+        pgi.reaction = "g6p_c --> f6p_c"
+        assert pgi.lower_bound == 0
+        pgi.bounds = (0, 1000)
+        assert pgi.bounds == (0, 1000)
+        assert not pgi.reversibility
+        pgi.reaction = "g6p_c <== f6p_c"
+        assert pgi.upper_bound == 0
+        assert pgi.reaction.strip() == "g6p_c <-- f6p_c"
+        pgi.reaction = "g6p_c --> f6p_c + h2o_c"
+        assert model.metabolites.h2o_c, pgi._metabolites
+        pgi.build_reaction_from_string("g6p_c --> f6p_c + foo", verbose=False)
+        assert model.metabolites.h2o_c not in pgi._metabolites
+        assert "foo" in model.metabolites
+        assert model.metabolites.foo in pgi._metabolites
+        assert len(model.metabolites) == m + 1
+
+    def test_copy(self, model):
+        PGI = model.reactions.PGI
+        copied = PGI.copy()
+        assert PGI is not copied
+        assert PGI._model is model
+        assert copied._model is not model
+        # the copy should refer to different metabolites and genes
+        for met in copied.metabolites:
+            assert met is not model.metabolites.get_by_id(met.id)
+            assert met.model is not model
+        for gene in copied.genes:
+            assert gene is not model.genes.get_by_id(gene.id)
+            assert gene.model is not model
+
+    def test_iadd(self, model):
+        PGI = model.reactions.PGI
+        EX_h2o = model.reactions.EX_h2o_e
+        original_PGI_gpr = PGI.gene_reaction_rule
+        PGI += EX_h2o
+        assert PGI.gene_reaction_rule == original_PGI_gpr
+        assert PGI.metabolites[model.metabolites.h2o_e] == -1.0
+        # original should not have changed
+        assert EX_h2o.gene_reaction_rule == ''
+        assert EX_h2o.metabolites[model.metabolites.h2o_e] == -1.0
+        # what about adding a reaction not in the model
+        new_reaction = Reaction("test")
+        new_reaction.add_metabolites({Metabolite("A"): -1, Metabolite("B"): 1})
+        PGI += new_reaction
+        assert PGI.gene_reaction_rule == original_PGI_gpr
+        assert len(PGI.gene_reaction_rule) == 5
+        # and vice versa
+        new_reaction += PGI
+        assert len(new_reaction.metabolites) == 5  # not
+        assert len(new_reaction.genes) == 1
+        assert new_reaction.gene_reaction_rule == original_PGI_gpr
+        # what about combining 2 gpr's
+        model.reactions.ACKr += model.reactions.ACONTa
+        expected_rule = '(b2296 or b3115 or b1849) and (b0118 or b1276)'
+        assert model.reactions.ACKr.gene_reaction_rule == expected_rule
+        assert len(model.reactions.ACKr.genes) == 5
+
+    def test_add(self, model):
+        # not in place addition should work on a copy
+        new = model.reactions.PGI + model.reactions.EX_h2o_e
+        assert new._model is not model
+        assert len(new.metabolites) == 3
+        # the copy should refer to different metabolites and genes
+        # This currently fails because add_metabolites does not copy.
+        # Should that be changed?
+        # for met in new.metabolites:
+        #    assert met is not model.metabolites.get_by_id(met.id)
+        #    assert met.model is not model
+        for gene in new.genes:
+            assert gene is not model.genes.get_by_id(gene.id)
+            assert gene.model is not model
+
+    def test_mul(self, model):
+        new = model.reactions.PGI * 2
+        assert set(new.metabolites.values()) == {-2, 2}
+
+    def test_sub(self, model):
+        new = model.reactions.PGI - model.reactions.EX_h2o_e
+        assert new._model is not model
+        assert len(new.metabolites) == 3
+
+
+class TestCobraMetabolites:
+    def test_metabolite_formula(self):
+        met = Metabolite("water")
+        met.formula = "H2O"
+        assert met.elements == {"H": 2, "O": 1}
+        assert met.formula_weight == 18.01528
+
+    def test_formula_element_setting(self, model):
+        met = model.metabolites[1]
+        orig_formula = str(met.formula)
+        orig_elements = dict(met.elements)
+        met.formula = ''
+        assert met.elements == {}
+        met.elements = orig_elements
+        assert met.formula == orig_formula
+
+
+class TestCobraModel:
+    """test core cobra functions"""
+
+    @pytest.mark.parametrize("solver", list(solver_dict))
+    def test_add_remove_reaction_benchmark(self, model, benchmark, solver):
+        metabolite_foo = Metabolite("test_foo")
+        metabolite_bar = Metabolite("test_bar")
+        metabolite_baz = Metabolite("test_baz")
+        actual_metabolite = model.metabolites[0]
+        dummy_reaction = Reaction("test_foo_reaction")
+        dummy_reaction.add_metabolites({metabolite_foo: -1,
+                                        metabolite_bar: 1,
+                                        metabolite_baz: -2,
+                                        actual_metabolite: 1})
+
+        def benchmark_add_reaction():
+            model.add_reaction(dummy_reaction)
+            if not getattr(model, 'solver', None):
+                solver_dict[solver].create_problem(model)
+            model.remove_reactions([dummy_reaction], delete=False)
+        benchmark(benchmark_add_reaction)
+
+    def test_add_reaction(self, model):
+        old_reaction_count = len(model.reactions)
+        old_metabolite_count = len(model.metabolites)
+        dummy_metabolite_1 = Metabolite("test_foo_1")
+        dummy_metabolite_2 = Metabolite("test_foo_2")
+        actual_metabolite = model.metabolites[0]
+        copy_metabolite = model.metabolites[1].copy()
+        dummy_reaction = Reaction("test_foo_reaction")
+        dummy_reaction.add_metabolites({dummy_metabolite_1: -1,
+                                        dummy_metabolite_2: 1,
+                                        copy_metabolite: -2,
+                                        actual_metabolite: 1})
+
+        model.add_reaction(dummy_reaction)
+        assert model.reactions.get_by_id(dummy_reaction.id) == dummy_reaction
+        for x in [dummy_metabolite_1, dummy_metabolite_2]:
+            assert model.metabolites.get_by_id(x.id) == x
+        # should have added 1 reaction and 2 metabolites
+        assert len(model.reactions) == old_reaction_count + 1
+        assert len(model.metabolites) == old_metabolite_count + 2
+        # tests on the added reaction
+        reaction_in_model = model.reactions.get_by_id(dummy_reaction.id)
+        assert type(reaction_in_model) is Reaction
+        assert reaction_in_model is dummy_reaction
+        assert len(reaction_in_model._metabolites) == 4
+        for i in reaction_in_model._metabolites:
+            assert type(i) == Metabolite
+        # tests on the added metabolites
+        met1_in_model = model.metabolites.get_by_id(dummy_metabolite_1.id)
+        assert met1_in_model is dummy_metabolite_1
+        copy_in_model = model.metabolites.get_by_id(copy_metabolite.id)
+        assert copy_metabolite is not copy_in_model
+        assert type(copy_in_model) is Metabolite
+        assert dummy_reaction in actual_metabolite._reaction
+        # test adding a different metabolite with the same name as an
+        # existing one uses the metabolite in the model
+        r2 = Reaction("test_foo_reaction2")
+        model.add_reaction(r2)
+        r2.add_metabolites({Metabolite(model.metabolites[0].id): 1})
+        assert model.metabolites[0] is list(r2._metabolites)[0]
+
+    def test_add_reaction_from_other_model(self, model):
+        other = model.copy()
+        for i in other.reactions:
+            i.id += "_other"
+        other.repair()
+        model.add_reactions(other.reactions)
+        # what if the other reaction has an error in its GPR
+        m1 = model.copy()
+        m2 = model.copy()
+        m1.reactions.PGI.remove_from_model()
+        m2.genes.b4025._reaction.clear()
+        m1.add_reaction(m2.reactions.PGI)
+
+    def test_model_remove_reaction(self, model):
+        old_reaction_count = len(model.reactions)
+        model.remove_reactions(["PGI"])
+        assert len(model.reactions) == old_reaction_count - 1
+        with pytest.raises(KeyError):
+            model.reactions.get_by_id("PGI")
+        model.remove_reactions(model.reactions[:1])
+        assert len(model.reactions) == old_reaction_count - 2
+        tmp_metabolite = Metabolite("testing")
+        model.reactions[0].add_metabolites({tmp_metabolite: 1})
+        assert tmp_metabolite in model.metabolites
+        model.remove_reactions(model.reactions[:1],
+                               remove_orphans=True)
+        assert tmp_metabolite not in model.metabolites
+
+    def test_reaction_remove(self, model):
+        old_reaction_count = len(model.reactions)
+        tmp_metabolite = Metabolite("testing")
+        # Delete without removing orphan
+        model.reactions[0].add_metabolites({tmp_metabolite: 1})
+        assert len(tmp_metabolite.reactions) == 1
+        # esnsure the stoichiometry is still the same using different objects
+        removed_reaction = model.reactions[0]
+        original_stoich = {i.id: value for i, value
+                           in removed_reaction._metabolites.items()}
+        model.reactions[0].remove_from_model(remove_orphans=False)
+        assert len(original_stoich) == len(removed_reaction._metabolites)
+        for met in removed_reaction._metabolites:
+            assert original_stoich[met.id] == removed_reaction._metabolites[
+                met]
+            assert met is not model.metabolites
+        # make sure it's still in the model
+        assert tmp_metabolite in model.metabolites
+        assert len(tmp_metabolite.reactions) == 0
+        assert len(model.reactions) == old_reaction_count - 1
+
+        # Now try it with removing orphans
+        model.reactions[0].add_metabolites({tmp_metabolite: 1})
+        assert len(tmp_metabolite.reactions) == 1
+        model.reactions[0].remove_from_model(remove_orphans=True)
+        assert tmp_metabolite not in model.metabolites
+        assert len(tmp_metabolite.reactions) == 0
+        assert len(model.reactions) == old_reaction_count - 2
+
+        # It shouldn't remove orphans if it's in 2 reactions however
+        model.reactions[0].add_metabolites({tmp_metabolite: 1})
+        model.reactions[1].add_metabolites({tmp_metabolite: 1})
+        assert len(tmp_metabolite.reactions) == 2
+        model.reactions[0].remove_from_model(remove_orphans=False)
+        assert tmp_metabolite in model.metabolites
+        assert len(tmp_metabolite.reactions) == 1
+        assert len(model.reactions) == old_reaction_count - 3
+
+    def test_reaction_delete(self, model):
+        old_reaction_count = len(model.reactions)
+        tmp_metabolite = Metabolite("testing")
+        # Delete without removing orphan
+        model.reactions[0].add_metabolites({tmp_metabolite: 1})
+        assert len(tmp_metabolite.reactions) == 1
+        model.reactions[0].delete(remove_orphans=False)
+        # make sure it's still in the model
+        assert tmp_metabolite in model.metabolites
+        assert len(tmp_metabolite.reactions) == 0
+        assert len(model.reactions) == old_reaction_count - 1
+
+        # Now try it with removing orphans
+        model.reactions[0].add_metabolites({tmp_metabolite: 1})
+        assert len(tmp_metabolite.reactions) == 1
+        model.reactions[0].delete(remove_orphans=True)
+        assert tmp_metabolite not in model.metabolites
+        assert len(tmp_metabolite.reactions) == 0
+        assert len(model.reactions) == old_reaction_count - 2
+
+        # It shouldn't remove orphans if it's in 2 reactions however
+        model.reactions[0].add_metabolites({tmp_metabolite: 1})
+        model.reactions[1].add_metabolites({tmp_metabolite: 1})
+        assert len(tmp_metabolite.reactions) == 2
+        model.reactions[0].delete(remove_orphans=False)
+        assert tmp_metabolite in model.metabolites
+        assert len(tmp_metabolite.reactions) == 1
+        assert len(model.reactions) == old_reaction_count - 3
+
+    def test_remove_gene(self, model):
+        target_gene = model.genes[0]
+        gene_reactions = list(target_gene.reactions)
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore")
+            target_gene.remove_from_model()
+        assert target_gene.model is None
+        # make sure the reaction was removed from the model
+        assert target_gene not in model.genes
+        # ensure the old reactions no longer have a record of the gene
+        for reaction in gene_reactions:
+            assert target_gene not in reaction.genes
+
+    @pytest.mark.parametrize("solver", list(solver_dict))
+    def test_copy_benchmark(self, model, solver, benchmark):
+        def _():
+            model.copy()
+            if not getattr(model, 'solver', None):
+                solver_dict[solver].create_problem(model)
+        benchmark(_)
+
+    @pytest.mark.parametrize("solver", list(solver_dict))
+    def test_copy_benchmark_large_model(self, large_model, solver, benchmark):
+        def _():
+            large_model.copy()
+            if not getattr(large_model, 'solver', None):
+                solver_dict[solver].create_problem(large_model)
+        benchmark(_)
+
+    def test_copy(self, model):
+        """modifying copy should not modify the original"""
+        # test that deleting reactions in the copy does not change the
+        # number of reactions in the original model
+        model_copy = model.copy()
+        old_reaction_count = len(model.reactions)
+        assert model_copy.notes is not model.notes
+        assert model_copy.annotation is not model.annotation
+        assert len(model.reactions) == len(model_copy.reactions)
+        assert len(model.metabolites) == len(model_copy.metabolites)
+        model_copy.remove_reactions(model_copy.reactions[0:5])
+        assert old_reaction_count == len(model.reactions)
+        assert len(model.reactions) != len(model_copy.reactions)
+
+    def test_deepcopy_benchmark(self, model, benchmark):
+        benchmark(deepcopy, model)
+
+    def test_deepcopy(self, model):
+        """Reference structures are maintained when deepcopying"""
+        model_copy = deepcopy(model)
+        for gene, gene_copy in zip(model.genes, model_copy.genes):
+            assert gene.id == gene_copy.id
+            reactions = sorted(i.id for i in gene.reactions)
+            reactions_copy = sorted(i.id for i in gene_copy.reactions)
+            assert reactions == reactions_copy
+        for reaction, reaction_copy in zip(model.reactions,
+                                           model_copy.reactions):
+            assert reaction.id == reaction_copy.id
+            metabolites = sorted(i.id for i in reaction._metabolites)
+            metabolites_copy = sorted(i.id for i in reaction_copy._metabolites)
+            assert metabolites == metabolites_copy
+
+    def test_add_reaction_orphans(self, model):
+        """test reaction addition
+
+        Need to verify that no orphan genes or metabolites are
+        contained in reactions after adding them to the model.
+        """
+        model = model.__class__('test')
+        model.add_reactions((x.copy() for x in model.reactions))
+        genes = []
+        metabolites = []
+        for x in model.reactions:
+            genes.extend(x.genes)
+            metabolites.extend(x._metabolites)
+        orphan_genes = [x for x in genes if x.model is not model]
+        orphan_metabolites = [x for x in metabolites if x.model is not model]
+        # check not dangling genes when running Model.add_reactions
+        assert len(orphan_genes) == 0
+        # 'check not dangling metabolites when running Model.add_reactions
+        assert len(orphan_metabolites) == 0
+
+    @pytest.mark.parametrize("solver", list(solver_dict))
+    def test_change_objective_benchmark(self, model, benchmark, solver):
+        atpm = model.reactions.get_by_id("ATPM")
+
+        def benchmark_change_objective():
+            model.objective = atpm.id
+            if not getattr(model, 'solver', None):
+                solver_dict[solver].create_problem(model)
+        benchmark(benchmark_change_objective)
+
+    def test_change_objective(self, model):
+        biomass = model.reactions.get_by_id("Biomass_Ecoli_core")
+        atpm = model.reactions.get_by_id("ATPM")
+        model.objective = atpm.id
+        assert atpm.objective_coefficient == 1.
+        assert biomass.objective_coefficient == 0.
+        assert model.objective == {atpm: 1.}
+        # change it back using object itself
+        model.objective = biomass
+        assert atpm.objective_coefficient == 0.
+        assert biomass.objective_coefficient == 1.
+        # set both to 1 with a list
+        model.objective = [atpm, biomass]
+        assert atpm.objective_coefficient == 1.
+        assert biomass.objective_coefficient == 1.
+        # set both using a dict
+        model.objective = {atpm: 0.2, biomass: 0.3}
+        assert atpm.objective_coefficient == 0.2
+        assert biomass.objective_coefficient == 0.3
+        # test setting by index
+        model.objective = model.reactions.index(atpm)
+        assert model.objective == {atpm: 1.}
+        # test by setting list of indexes
+        model.objective = map(model.reactions.index, [atpm, biomass])
+        assert model.objective == {atpm: 1., biomass: 1.}
+
+
+ at pytest.mark.skipif(scipy is None, reason="scipy required for ArrayBasedModel")
+class TestCobraArrayModel:
+    def test_array_model(self, model):
+        for matrix_type in ["scipy.dok_matrix", "scipy.lil_matrix"]:
+            array_model = model.to_array_based_model(matrix_type=matrix_type)
+            assert array_model.S[7, 0] == -1
+            assert array_model.S[43, 0] == 0
+            array_model.S[43, 0] = 1
+            assert array_model.S[43, 0] == 1
+            assert array_model.reactions[0]._metabolites[
+                       array_model.metabolites[43]] == 1
+            array_model.S[43, 0] = 0
+            assert array_model.lower_bounds[0] == array_model.reactions[
+                0].lower_bound
+            assert array_model.lower_bounds[5] == array_model.reactions[
+                5].lower_bound
+            assert array_model.upper_bounds[0] == array_model.reactions[
+                0].upper_bound
+            assert array_model.upper_bounds[5] == array_model.reactions[
+                5].upper_bound
+            array_model.lower_bounds[6] = 2
+            assert array_model.lower_bounds[6] == 2
+            assert array_model.reactions[6].lower_bound == 2
+            # this should fail because it is the wrong size
+            with pytest.raises(Exception):
+                array_model.upper_bounds = [0, 1]
+            array_model.upper_bounds = [0] * len(array_model.reactions)
+            assert max(array_model.upper_bounds) == 0
+            # test something for all the attributes
+            array_model.lower_bounds[2] = -1
+            assert array_model.reactions[2].lower_bound == -1
+            assert array_model.lower_bounds[2] == -1
+            array_model.objective_coefficients[2] = 1
+            assert array_model.reactions[2].objective_coefficient == 1
+            assert array_model.objective_coefficients[2] == 1
+            array_model.b[2] = 1
+            assert array_model.metabolites[2]._bound == 1
+            assert array_model.b[2] == 1
+            array_model.constraint_sense[2] = "L"
+            assert array_model.metabolites[2]._constraint_sense == "L"
+            assert array_model.constraint_sense[2] == "L"
+            # test resize matrix on reaction removal
+            m, n = array_model.S.shape
+            array_model.remove_reactions([array_model.reactions[2]],
+                                         remove_orphans=False)
+            assert len(array_model.metabolites) == array_model.S.shape[0]
+            assert len(array_model.reactions) == array_model.S.shape[1]
+            assert array_model.S.shape == (m, n - 1)
+
+    def test_array_based_model_add(self, model):
+        array_model = model.to_array_based_model()
+        m = len(array_model.metabolites)
+        n = len(array_model.reactions)
+        for matrix_type in ["scipy.dok_matrix", "scipy.lil_matrix"]:
+            test_model = model.copy().to_array_based_model(
+                matrix_type=matrix_type)
+            test_reaction = Reaction("test")
+            test_reaction.add_metabolites({test_model.metabolites[0]: 4})
+            test_reaction.lower_bound = -3.14
+            test_model.add_reaction(test_reaction)
+            assert len(test_model.reactions) == n + 1
+            assert test_model.S.shape == (m, n + 1)
+            assert len(test_model.lower_bounds) == n + 1
+            assert len(test_model.upper_bounds) == n + 1
+            assert test_model.S[0, n] == 4
+            assert test_model.S[7, 0] == -1
+            assert test_model.lower_bounds[n] == -3.14
+
+    def test_array_based_select(self, array_model):
+        atpm_select = array_model.reactions[array_model.lower_bounds > 0]
+        assert len(atpm_select) == 1
+        assert atpm_select[0].id == "ATPM"
+        assert len(
+            array_model.reactions[array_model.lower_bounds <= 0]) == len(
+            array_model.reactions) - 1
+        # mismatched dimensions should give an error
+        with pytest.raises(TypeError):
+            array_model.reactions[[True, False]]
+
+    def test_array_based_bounds_setting(self, array_model):
+        model = array_model
+        bounds = [0.0] * len(model.reactions)
+        model.lower_bounds = bounds
+        assert type(model.reactions[0].lower_bound) == float
+        assert abs(model.reactions[0].lower_bound) < 10 ** -5
+        model.upper_bounds[1] = 1234.0
+        assert abs(model.reactions[1].upper_bound - 1234.0) < 10 ** -5
+        model.upper_bounds[9:11] = [100.0, 200.0]
+        assert abs(model.reactions[9].upper_bound - 100.0) < 10 ** -5
+        assert abs(model.reactions[10].upper_bound - 200.0) < 10 ** -5
+        model.upper_bounds[9:11] = 123.0
+        assert abs(model.reactions[9].upper_bound - 123.0) < 10 ** -5
+        assert abs(model.reactions[10].upper_bound - 123.0) < 10 ** -5
diff --git a/cobra/test/test_solvers.py b/cobra/test/test_solvers.py
new file mode 100644
index 0000000..4b32eab
--- /dev/null
+++ b/cobra/test/test_solvers.py
@@ -0,0 +1,266 @@
+import pytest
+from cobra.core import Model, Reaction, Metabolite
+from cobra import solvers
+from .conftest import model
+
+try:
+    import scipy
+except ImportError:
+    scipy = None
+
+
+ at pytest.fixture(scope="class", params=list(solvers.solver_dict))
+def solver_test(request):
+    solver = solvers.solver_dict[request.param]
+    old_solution = 0.8739215
+    infeasible_model = Model()
+    metabolite_1 = Metabolite("met1")
+    reaction_1 = Reaction("rxn1")
+    reaction_2 = Reaction("rxn2")
+    reaction_1.add_metabolites({metabolite_1: 1})
+    reaction_2.add_metabolites({metabolite_1: 1})
+    reaction_1.lower_bound = 1
+    reaction_2.upper_bound = 2
+    infeasible_model.add_reactions([reaction_1, reaction_2])
+    return solver, old_solution, infeasible_model
+
+
+class TestCobraSolver:
+    def test_attributes(self, solver_test):
+        solver, old_solution, infeasible_model = solver_test
+        assert hasattr(solver, "create_problem")
+        assert hasattr(solver, "solve_problem")
+        assert hasattr(solver, "get_status")
+        assert hasattr(solver, "get_objective_value")
+        assert hasattr(solver, "format_solution")
+        assert hasattr(solver, "change_variable_bounds")
+        assert hasattr(solver, "change_variable_objective")
+        assert hasattr(solver, "solve")
+        assert hasattr(solver, "set_parameter")
+
+    def test_creation(self, solver_test, model):
+        solver, old_solution, infeasible_model = solver_test
+        solver.create_problem(model)
+
+    def test_solve_feasible(self, solver_test, model):
+        solver, old_solution, infeasible_model = solver_test
+        solution = solver.solve(model)
+        assert solution.status == "optimal"
+        assert abs(old_solution - solution.f) < 10 ** -4
+
+    def test_solve_minimize(self, solver_test, model):
+        solver, old_solution, infeasible_model = solver_test
+        solution = solver.solve(model, objective_sense='minimize')
+        assert solution.status == "optimal"
+        assert abs(solution.f) < 10 ** -4
+
+    def test_low_level_control(self, solver_test):
+        solver, old_solution, infeasible_model = solver_test
+        lp = solver.create_problem(infeasible_model)
+        solver.solve_problem(lp)
+        assert solver.get_status(lp) == "infeasible"
+        # going to make feasible
+        solver.change_variable_bounds(lp, 0, -2., 2.)
+        solver.change_variable_bounds(lp, 1, -2., 2.)
+        solver.solve_problem(lp)
+        # should now be feasible, but obj = 0
+        assert solver.get_status(lp) == "optimal"
+        assert abs(solver.get_objective_value(lp)) < 10 ** -4
+        # should now have obj = 2 (maximize should be the default)
+        solver.change_variable_objective(lp, 0, 1.)
+        solver.solve_problem(lp)
+        assert abs(solver.get_objective_value(lp) - 2) < 10 ** -4
+        # should now solve with obj = -2
+        solver.solve_problem(lp, objective_sense="minimize")
+        assert abs(solver.get_objective_value(lp) + 2) < 10 ** -4
+        # should now have obj = 4
+        solver.change_variable_objective(lp, 0, 2.)
+        solver.solve_problem(lp, objective_sense="maximize")
+        assert abs(solver.get_objective_value(lp) - 4) < 10 ** -4
+        # make sure the solution looks good still
+        solution = solver.format_solution(lp, infeasible_model)
+        assert abs(solution.x[0] - 2) < 10 ** -4
+        assert abs(solution.x[1] + 2) < 10 ** -4
+        assert abs(solution.x_dict["rxn1"] - 2) < 10 ** -4
+        assert abs(solution.x_dict["rxn2"] + 2) < 10 ** -4
+
+    def test_set_objective_sense(self, solver_test, model):
+        solver, old_solution, infeasible_model = solver_test
+        maximize = solver.create_problem(model, objective_sense="maximize")
+        minimize = solver.create_problem(model, objective_sense="minimize")
+        solver.solve_problem(maximize)
+        solver.solve_problem(minimize)
+        max_solution = solver.format_solution(maximize, model)
+        min_solution = solver.format_solution(minimize, model)
+        assert min_solution.status == "optimal"
+        assert abs(min_solution.f) < 10 ** -4
+        assert abs(old_solution - max_solution.f) < 10 ** -4
+        assert max_solution.status == "optimal"
+        # if we set minimize at creation, can we override it at solve
+        solver.solve_problem(minimize, objective_sense="maximize")
+        override_minimize = solver.format_solution(minimize, model)
+        assert abs(max_solution.f - override_minimize.f) < 10 ** -4
+
+    def test_solve_mip(self, solver_test):
+        solver, old_solution, infeasible_model = solver_test
+        if not hasattr(solver, "_SUPPORTS_MILP") or not solver._SUPPORTS_MILP:
+            pytest.skip("no milp support")
+        cobra_model = Model('MILP_implementation_test')
+        constraint = Metabolite("constraint")
+        constraint._bound = 2.5
+        x = Reaction("x")
+        x.lower_bound = 0.
+        x.objective_coefficient = 1.
+        x.add_metabolites({constraint: 2.5})
+        y = Reaction("y")
+        y.lower_bound = 0.
+        y.objective_coefficient = 1.
+        y.add_metabolites({constraint: 1.})
+        cobra_model.add_reactions([x, y])
+        float_sol = solver.solve(cobra_model)
+        # add an integer constraint
+        y.variable_kind = "integer"
+        int_sol = solver.solve(cobra_model)
+        assert abs(float_sol.f - 2.5) < 10 ** -5
+        assert abs(float_sol.x_dict["y"] - 2.5) < 10 ** -5
+        assert int_sol.status == "optimal"
+        assert abs(int_sol.f - 2.2) < 10 ** -3
+        assert abs(int_sol.x_dict["y"] - 2.0) < 10 ** -3
+
+    def test_solve_infeasible(self, solver_test):
+        solver, old_solution, infeasible_model = solver_test
+        solution = solver.solve(infeasible_model)
+        assert solution.status == "infeasible"
+
+    def test_independent_creation(self, solver_test, model):
+        solver, old_solution, infeasible_model = solver_test
+        feasible_lp = solver.create_problem(model)
+        infeasible_lp = solver.create_problem(infeasible_model)
+        solver.solve_problem(feasible_lp)
+        solver.solve_problem(infeasible_lp)
+        feasible_solution = solver.format_solution(feasible_lp, model)
+        infeasible_solution = solver.format_solution(infeasible_lp,
+                                                     infeasible_model)
+        assert feasible_solution.status == "optimal"
+        assert abs(old_solution - feasible_solution.f) < 10 ** -4
+        assert infeasible_solution.status == "infeasible"
+
+    def test_change_coefficient(self, solver_test):
+        solver, old_solution, infeasible_model = solver_test
+        c = Metabolite("c")
+        c._bound = 6
+        x = Reaction("x")
+        x.lower_bound = 1.
+        y = Reaction("y")
+        y.lower_bound = 0.
+        x.add_metabolites({c: 1})
+        z = Reaction("z")
+        z.add_metabolites({c: 1})
+        z.objective_coefficient = 1
+        m = Model("test_model")
+        m.add_reactions([x, y, z])
+        # change an existing coefficient
+        lp = solver.create_problem(m)
+        solver.solve_problem(lp)
+        sol1 = solver.format_solution(lp, m)
+        assert sol1.status == "optimal"
+        solver.change_coefficient(lp, 0, 0, 2)
+        solver.solve_problem(lp)
+        sol2 = solver.format_solution(lp, m)
+        assert sol2.status == "optimal"
+        assert abs(sol1.f - 5.0) < 10 ** -3
+        assert abs(sol2.f - 4.0) < 10 ** -3
+        # change a new coefficient
+        z.objective_coefficient = 0.
+        y.objective_coefficient = 1.
+        lp = solver.create_problem(m)
+        solver.change_coefficient(lp, 0, 1, 2)
+        solver.solve_problem(lp)
+        solution = solver.format_solution(lp, m)
+        assert solution.status == "optimal"
+        assert abs(solution.x_dict["y"] - 2.5) < 10 ** -3
+
+    def test_inequality(self, solver_test):
+        solver, old_solution, infeasible_model = solver_test
+        # The space enclosed by the constraints is a 2D triangle with
+        # vertexes as (3, 0), (1, 2), and (0, 1)
+        # c1 encodes y - x > 1 ==> y > x - 1
+        # c2 encodes y + x < 3 ==> y < 3 - x
+        c1 = Metabolite("c1")
+        c2 = Metabolite("c2")
+        x = Reaction("x")
+        x.lower_bound = 0
+        y = Reaction("y")
+        y.lower_bound = 0
+        x.add_metabolites({c1: -1, c2: 1})
+        y.add_metabolites({c1: 1, c2: 1})
+        c1._bound = 1
+        c1._constraint_sense = "G"
+        c2._bound = 3
+        c2._constraint_sense = "L"
+        m = Model()
+        m.add_reactions([x, y])
+        # test that optimal values are at the vertices
+        m.objective = "x"
+        assert abs(solver.solve(m).f - 1.0) < 10 ** -3
+        assert abs(solver.solve(m).x_dict["y"] - 2.0) < 10 ** -3
+        m.objective = "y"
+        assert abs(solver.solve(m).f - 3.0) < 10 ** -3
+        assert abs(
+            solver.solve(m, objective_sense="minimize").f - 1.0) < 10 ** -3
+
+    @pytest.mark.skipif(scipy is None,
+                        reason="scipy required for quadratic objectives")
+    def test_quadratic(self, solver_test):
+        solver, old_solution, infeasible_model = solver_test
+        if not hasattr(solver, "set_quadratic_objective"):
+            pytest.skip("no qp support")
+        c = Metabolite("c")
+        c._bound = 2
+        x = Reaction("x")
+        x.objective_coefficient = -0.5
+        x.lower_bound = 0.
+        y = Reaction("y")
+        y.objective_coefficient = -0.5
+        y.lower_bound = 0.
+        x.add_metabolites({c: 1})
+        y.add_metabolites({c: 1})
+        m = Model()
+        m.add_reactions([x, y])
+        lp = solver.create_problem(m)
+        quadratic_obj = scipy.sparse.eye(2) * 2
+        solver.set_quadratic_objective(lp, quadratic_obj)
+        solver.solve_problem(lp, objective_sense="minimize")
+        solution = solver.format_solution(lp, m)
+        assert solution.status == "optimal"
+        # Respecting linear objectives also makes the objective value 1.
+        assert abs(solution.f - 1.) < 10 ** -3
+        assert abs(solution.x_dict["y"] - 1.) < 10 ** -3
+        assert abs(solution.x_dict["y"] - 1.) < 10 ** -3
+        # When the linear objectives are removed the objective value is 2.
+        solver.change_variable_objective(lp, 0, 0.)
+        solver.change_variable_objective(lp, 1, 0.)
+        solver.solve_problem(lp, objective_sense="minimize")
+        solution = solver.format_solution(lp, m)
+        assert solution.status == "optimal"
+        assert abs(solution.f - 2.) < 10 ** -3
+        # test quadratic from solve function
+        solution = solver.solve(m, quadratic_component=quadratic_obj,
+                                objective_sense="minimize")
+        assert solution.status == "optimal"
+        assert abs(solution.f - 1.) < 10 ** -3
+        c._bound = 6
+        z = Reaction("z")
+        x.objective_coefficient = 0.
+        y.objective_coefficient = 0.
+        z.lower_bound = 0.
+        z.add_metabolites({c: 1})
+        m.add_reaction(z)
+        solution = solver.solve(m, quadratic_component=scipy.sparse.eye(3),
+                                objective_sense="minimize")
+        # should be 12 not 24 because 1/2 (V^T Q V)
+        assert solution.status == "optimal"
+        assert abs(solution.f - 6) < 10 ** -3
+        assert abs(solution.x_dict["x"] - 2) < 10 ** -6
+        assert abs(solution.x_dict["y"] - 2) < 10 ** -6
+        assert abs(solution.x_dict["z"] - 2) < 10 ** -6
diff --git a/cobra/test/unit_tests.py b/cobra/test/unit_tests.py
deleted file mode 100644
index fdd583b..0000000
--- a/cobra/test/unit_tests.py
+++ /dev/null
@@ -1,806 +0,0 @@
-import sys
-from unittest import TestCase, TestLoader, TextTestRunner, skipIf
-from copy import copy, deepcopy
-from pickle import loads, dumps, HIGHEST_PROTOCOL
-import warnings
-import re
-
-if __name__ == "__main__":
-    sys.path.insert(0, "../..")
-    from cobra.test import create_test_model
-    from cobra import Object, Model, Metabolite, Reaction, DictList
-    sys.path.pop(0)
-else:
-    from . import create_test_model
-    from .. import Object, Model, Metabolite, Reaction, DictList
-
-# libraries which may or may not be installed
-libraries = ["scipy"]
-for library in libraries:
-    try:
-        exec("import %s" % library)
-    except ImportError:
-        exec("%s = None" % library)
-
-
-class TestDictList(TestCase):
-    def setUp(self):
-        self.obj = Object("test1")
-        self.list = DictList()
-        self.list.append(self.obj)
-
-    def testContains(self):
-        self.assertIn(self.obj, self.list)
-        self.assertIn(self.obj.id, self.list)
-        self.assertNotIn(Object("not_in"), self.list)
-        self.assertNotIn('not_in', self.list)
-
-    def testIndex(self):
-        self.assertEqual(self.list.index("test1"), 0)
-        self.assertEqual(self.list.index(self.obj), 0)
-        self.assertRaises(ValueError, self.list.index, "f")
-        self.assertRaises(ValueError, self.list.index, Object("f"))
-        # ensure trying to index with an object that is a different object
-        # also raises an error
-        self.assertRaises(ValueError, self.list.index, Object("test1"))
-
-    def testIndependent(self):
-        a = DictList([Object("o1"), Object("o2")])
-        b = DictList()
-        self.assertIn("o1", a)
-        self.assertNotIn("o1", b)
-        b.append(Object("o3"))
-        self.assertNotIn("o3", a)
-        self.assertIn("o3", b)
-
-    def testAppend(self):
-        obj2 = Object("test2")
-        self.list.append(obj2)
-        self.assertRaises(ValueError, self.list.append, Object("test1"))
-        self.assertEqual(self.list.index(obj2), 1)
-        self.assertEqual(self.list[1], obj2)
-        self.assertIs(self.list.get_by_id("test2"), obj2)
-        self.assertEqual(len(self.list), 2)
-
-    def testInsert(self):
-        obj2 = Object("a")
-        self.list.insert(0, obj2)
-        self.assertEqual(self.list.index(obj2), 0)
-        self.assertEqual(self.list.index("test1"), 1)
-        self.assertIs(self.list.get_by_id("a"), obj2)
-        self.assertEqual(len(self.list), 2)
-        self.assertRaises(ValueError, self.list.append, obj2)
-
-    def testExtend(self):
-        obj_list = [Object("test%d" % (i)) for i in range(2, 10)]
-        self.list.extend(obj_list)
-        self.assertEqual(self.list[1].id, "test2")
-        self.assertEqual(self.list.get_by_id("test2"), obj_list[0])
-        self.assertEqual(self.list[8].id, "test9")
-        self.assertEqual(len(self.list), 9)
-        self.assertRaises(ValueError, self.list.extend, [Object("test1")])
-        # Even if the object is unique, if it is present twice in the new
-        # list, it should still raise an exception
-        self.assertRaises(ValueError, self.list.extend,
-                          [Object("testd"), Object("testd")])
-
-    def testIadd(self):
-        obj_list = [Object("test%d" % (i)) for i in range(2, 10)]
-        self.list += obj_list
-        self.assertEqual(self.list[1].id, "test2")
-        self.assertEqual(self.list.get_by_id("test2"), obj_list[0])
-        self.assertEqual(self.list[8].id, "test9")
-        self.assertEqual(len(self.list), 9)
-
-    def testAdd(self):
-        obj_list = [Object("test%d" % (i)) for i in range(2, 10)]
-        sum = self.list + obj_list
-        self.assertIsNot(sum, self.list)
-        self.assertIsNot(sum, obj_list)
-        self.assertEqual(self.list[0].id, "test1")
-        self.assertEqual(sum[1].id, "test2")
-        self.assertEqual(sum.get_by_id("test2"), obj_list[0])
-        self.assertEqual(sum[8].id, "test9")
-        self.assertEqual(len(self.list), 1)
-        self.assertEqual(len(sum), 9)
-
-    def testInitCopy(self):
-        self.list.append(Object("test2"))
-        copied = DictList(self.list)
-        self.assertIsNot(self.list, copied)
-        self.assertIsInstance(copied, self.list.__class__)
-        self.assertEqual(len(self.list), len(copied))
-        for i, v in enumerate(self.list):
-            self.assertEqual(self.list[i].id, copied[i].id)
-            self.assertEqual(i, copied.index(v.id))
-            self.assertIs(self.list[i], copied[i])
-            self.assertIs(v, copied.get_by_id(v.id))
-
-    def testSlice(self):
-        self.list.append(Object("test2"))
-        self.list.append(Object("test3"))
-        sliced = self.list[:-1]
-        self.assertIsNot(self.list, sliced)
-        self.assertIsInstance(sliced, self.list.__class__)
-        self.assertEqual(len(self.list), len(sliced) + 1)
-        for i, v in enumerate(sliced):
-            self.assertEqual(self.list[i].id, sliced[i].id)
-            self.assertEqual(i, sliced.index(v.id))
-            self.assertIs(self.list[i], sliced[i])
-            self.assertIs(self.list[i], sliced.get_by_id(v.id))
-
-    def testCopy(self):
-        self.list.append(Object("test2"))
-        copied = copy(self.list)
-        self.assertIsNot(self.list, copied)
-        self.assertIsInstance(copied, self.list.__class__)
-        self.assertEqual(len(self.list), len(copied))
-        for i, v in enumerate(self.list):
-            self.assertEqual(self.list[i].id, copied[i].id)
-            self.assertEqual(i, copied.index(v.id))
-            self.assertIs(self.list[i], copied[i])
-            self.assertIs(v, copied.get_by_id(v.id))
-
-    def testDeepcopy(self):
-        self.list.append(Object("test2"))
-        copied = deepcopy(self.list)
-        self.assertIsNot(self.list, copied)
-        self.assertIsInstance(copied, self.list.__class__)
-        self.assertEqual(len(self.list), len(copied))
-        for i, v in enumerate(self.list):
-            self.assertEqual(self.list[i].id, copied[i].id)
-            self.assertEqual(i, copied.index(v.id))
-            self.assertIsNot(self.list[i], copied[i])
-            self.assertIsNot(v, copied.get_by_id(v.id))
-
-    def testPickle(self):
-        self.list.append(Object("test2"))
-        for protocol in range(HIGHEST_PROTOCOL):
-            pickle_str = dumps(self.list, protocol=protocol)
-            copied = loads(pickle_str)
-            self.assertIsNot(self.list, copied)
-            self.assertIsInstance(copied, self.list.__class__)
-            self.assertEqual(len(self.list), len(copied))
-            for i, v in enumerate(self.list):
-                self.assertEqual(self.list[i].id, copied[i].id)
-                self.assertEqual(i, copied.index(v.id))
-                self.assertIsNot(self.list[i], copied[i])
-                self.assertIsNot(v, copied.get_by_id(v.id))
-
-    def testQuery(self):
-        obj2 = Object("test2")
-        obj2.name = "foobar1"
-        self.list.append(obj2)
-        result = self.list.query("test1")  # matches only test1
-        self.assertEqual(len(result), 1)
-        self.assertEqual(result[0], self.obj)
-        result = self.list.query("foo", "name")  # matches only test2
-        self.assertEqual(len(result), 1)
-        self.assertEqual(result[0], obj2)
-        result = self.list.query("test")  # matches test1 and test2
-        self.assertEqual(len(result), 2)
-        # test with a regular expression
-        result = self.list.query(re.compile("test[0-9]"))
-        self.assertEqual(len(result), 2)
-        result = self.list.query(re.compile("test[29]"))
-        self.assertEqual(len(result), 1)
-        # test query of name
-        result = self.list.query(re.compile("foobar."), "name")
-        self.assertEqual(len(result), 1)
-
-    def testRemoval(self):
-        obj_list = DictList(Object("test%d" % (i)) for i in range(2, 10))
-        del obj_list[3]
-        self.assertNotIn("test5", obj_list)
-        self.assertEqual(obj_list.index(obj_list[-1]), len(obj_list) - 1)
-        self.assertEqual(len(obj_list), 7)
-        del obj_list[3:5]
-        self.assertNotIn("test6", obj_list)
-        self.assertNotIn("test7", obj_list)
-        self.assertEqual(obj_list.index(obj_list[-1]), len(obj_list) - 1)
-        self.assertEqual(len(obj_list), 5)
-        removed = obj_list.pop(1)
-        self.assertEqual(obj_list.index(obj_list[-1]), len(obj_list) - 1)
-        self.assertEqual(removed.id, "test3")
-        self.assertNotIn("test3", obj_list)
-        self.assertEqual(len(obj_list), 4)
-        removed = obj_list.pop()
-        self.assertEqual(removed.id, "test9")
-        self.assertNotIn(removed.id, obj_list)
-        self.assertEqual(len(obj_list), 3)
-
-    def testSet(self):
-        obj_list = DictList(Object("test%d" % (i)) for i in range(10))
-        obj_list[4] = Object("testa")
-        self.assertEqual(obj_list.index("testa"), 4)
-        self.assertEqual(obj_list[4].id, "testa")
-        obj_list[5:7] = [Object("testb"), Object("testc")]
-        self.assertEqual(obj_list.index("testb"), 5)
-        self.assertEqual(obj_list[5].id, "testb")
-        self.assertEqual(obj_list.index("testc"), 6)
-        self.assertEqual(obj_list[6].id, "testc")
-        # Even if the object is unique, if it is present twice in the new
-        # list, it should still raise an exception
-        self.assertRaises(ValueError, obj_list.__setitem__, slice(5, 7),
-                          [Object("testd"), Object("testd")])
-
-    def testSortandReverse(self):
-        dl = DictList(Object("test%d" % (i)) for i in reversed(range(10)))
-        self.assertEqual(dl[0].id, "test9")
-        dl.sort()
-        self.assertEqual(len(dl), 10)
-        self.assertEqual(dl[0].id, "test0")
-        self.assertEqual(dl.index("test0"), 0)
-        dl.reverse()
-        self.assertEqual(dl[0].id, "test9")
-        self.assertEqual(dl.index("test0"), 9)
-
-    def testDir(self):
-        """makes sure tab complete will work"""
-        attrs = dir(self.list)
-        self.assertIn("test1", attrs)
-        self.assertIn("_dict", attrs)  # attribute of DictList
-
-    def testUnion(self):
-        self.list.union([Object("test1"), Object("test2")])
-        # should only add 1 element
-        self.assertEqual(len(self.list), 2)
-        self.assertEqual(self.list.index("test2"), 1)
-
-
-class CobraTestCase(TestCase):
-    def setUp(self):
-        self.model = create_test_model("textbook")
-        self.model_class = Model
-
-
-class TestReactions(CobraTestCase):
-    def testGPR(self):
-        model = self.model_class()
-        reaction = Reaction("test")
-        # set a gpr to  reaction not in a model
-        reaction.gene_reaction_rule = "(g1 or g2) and g3"
-        self.assertEqual(reaction.gene_reaction_rule, "(g1 or g2) and g3")
-        self.assertEqual(len(reaction.genes), 3)
-        # adding reaction with a GPR propagates to the model
-        model.add_reaction(reaction)
-        self.assertEqual(len(model.genes), 3)
-        # ensure the gene objects are the same in the model and reaction
-        reaction_gene = list(reaction.genes)[0]
-        model_gene = model.genes.get_by_id(reaction_gene.id)
-        self.assertIs(reaction_gene, model_gene)
-        # test ability to handle uppercase AND/OR
-        with warnings.catch_warnings():
-            warnings.simplefilter("ignore")
-            reaction.gene_reaction_rule = "(b1 AND b2) OR (b3 and b4)"
-        self.assertEqual(reaction.gene_reaction_rule,
-                         "(b1 and b2) or (b3 and b4)")
-        self.assertEqual(len(reaction.genes), 4)
-        # ensure regular expressions correctly extract genes from malformed
-        # GPR string
-        with warnings.catch_warnings():
-            warnings.simplefilter("ignore")
-            reaction.gene_reaction_rule = "(a1 or a2"
-            self.assertEqual(len(reaction.genes), 2)
-            reaction.gene_reaction_rule = "(forT or "
-            self.assertEqual(len(reaction.genes), 1)
-
-    def testGPR_modification(self):
-        model = self.model
-        reaction = model.reactions.get_by_id("PGI")
-        old_gene = list(reaction.genes)[0]
-        new_gene = model.genes.get_by_id("s0001")
-        # add an existing 'gene' to the gpr
-        reaction.gene_reaction_rule = 's0001'
-        self.assertIn(new_gene, reaction.genes)
-        self.assertIn(reaction, new_gene.reactions)
-        # removed old gene correctly
-        self.assertNotIn(old_gene, reaction.genes)
-        self.assertNotIn(reaction, old_gene.reactions)
-        # add a new 'gene' to the gpr
-        reaction.gene_reaction_rule = 'fake_gene'
-        self.assertTrue(model.genes.has_id("fake_gene"))
-        fake_gene = model.genes.get_by_id("fake_gene")
-        self.assertIn(fake_gene, reaction.genes)
-        self.assertIn(reaction, fake_gene.reactions)
-        fake_gene.name = "foo_gene"
-        self.assertEqual(reaction.gene_name_reaction_rule, fake_gene.name)
-
-    def test_add_metabolite(self):
-        model = self.model
-        reaction = model.reactions.get_by_id("PGI")
-        reaction.add_metabolites({model.metabolites[0]: 1})
-        self.assertIn(model.metabolites[0], reaction._metabolites)
-        fake_metabolite = Metabolite("fake")
-        reaction.add_metabolites({fake_metabolite: 1})
-        self.assertIn(fake_metabolite, reaction._metabolites)
-        self.assertTrue(model.metabolites.has_id("fake"))
-        self.assertIs(model.metabolites.get_by_id("fake"), fake_metabolite)
-
-        # test adding by string
-        reaction.add_metabolites({"g6p_c": -1})  # already in reaction
-        self.assertTrue(
-            reaction._metabolites[model.metabolites.get_by_id("g6p_c")], -2)
-        reaction.add_metabolites({"h_c": 1})  # not currently in reaction
-        self.assertTrue(
-            reaction._metabolites[model.metabolites.get_by_id("h_c")], 1)
-        with self.assertRaises(KeyError):
-            reaction.add_metabolites({"missing": 1})
-
-        # test adding to a new Reaction
-        reaction = Reaction("test")
-        self.assertEqual(len(reaction._metabolites), 0)
-        reaction.add_metabolites({Metabolite("test_met"): -1})
-        self.assertEqual(len(reaction._metabolites), 1)
-
-    def test_subtract_metabolite(self):
-        model = self.model
-        reaction = model.reactions.get_by_id("PGI")
-        reaction.subtract_metabolites(reaction.metabolites)
-        self.assertEqual(len(reaction.metabolites), 0)
-
-    def test_mass_balance(self):
-        model = self.model
-        reaction = model.reactions.get_by_id("PGI")
-        # should be balanced now
-        self.assertEqual(len(reaction.check_mass_balance()), 0)
-        # should not be balanced after adding a hydrogen
-        reaction.add_metabolites({model.metabolites.get_by_id("h_c"): 1})
-        imbalance = reaction.check_mass_balance()
-        self.assertIn("H", imbalance)
-        self.assertEqual(imbalance["H"], 1)
-
-    def test_build_from_string(self):
-        model = self.model
-        m = len(model.metabolites)
-        pgi = model.reactions.get_by_id("PGI")
-        pgi.reaction = "g6p_c --> f6p_c"
-        self.assertEqual(pgi.lower_bound, 0)
-        pgi.bounds = (0, 1000)
-        self.assertEqual(pgi.bounds, (0, 1000))
-        self.assertEqual(pgi.reversibility, False)
-        pgi.reaction = "g6p_c <== f6p_c"
-        self.assertEqual(pgi.upper_bound, 0)
-        self.assertEqual(pgi.reaction.strip(), "g6p_c <-- f6p_c")
-        pgi.reaction = "g6p_c --> f6p_c + h2o_c"
-        self.assertIn(model.metabolites.h2o_c, pgi._metabolites)
-        pgi.build_reaction_from_string("g6p_c --> f6p_c + foo", verbose=False)
-        self.assertNotIn(model.metabolites.h2o_c, pgi._metabolites)
-        self.assertIn("foo", model.metabolites)
-        self.assertIn(model.metabolites.foo, pgi._metabolites)
-        self.assertEqual(len(model.metabolites), m + 1)
-
-    def test_copy(self):
-        model = self.model
-        PGI = model.reactions.PGI
-        copied = PGI.copy()
-        self.assertIsNot(PGI, copied)
-        self.assertIs(PGI._model, model)
-        self.assertIsNot(copied._model, model)
-        # the copy should refer to different metabolites and genes
-        for met in copied.metabolites:
-            self.assertIsNot(met, model.metabolites.get_by_id(met.id))
-            self.assertIsNot(met.model, model)
-        for gene in copied.genes:
-            self.assertIsNot(gene, model.genes.get_by_id(gene.id))
-            self.assertIsNot(gene.model, model)
-
-    def test_iadd(self):
-        model = self.model
-        PGI = model.reactions.PGI
-        EX_h2o = model.reactions.EX_h2o_e
-        original_PGI_gpr = PGI.gene_reaction_rule
-        PGI += EX_h2o
-        self.assertEqual(PGI.gene_reaction_rule, original_PGI_gpr)
-        self.assertEqual(PGI.metabolites[model.metabolites.h2o_e], -1.0)
-        # original should not have changed
-        self.assertEqual(EX_h2o.gene_reaction_rule, '')
-        self.assertEqual(EX_h2o.metabolites[model.metabolites.h2o_e], -1.0)
-        # what about adding a reaction not in the model
-        new_reaction = Reaction("test")
-        new_reaction.add_metabolites({Metabolite("A"): -1, Metabolite("B"): 1})
-        PGI += new_reaction
-        self.assertEqual(PGI.gene_reaction_rule, original_PGI_gpr)
-        self.assertEqual(len(PGI.gene_reaction_rule), 5)
-        # and vice versa
-        new_reaction += PGI
-        self.assertEqual(len(new_reaction.metabolites), 5)  # not 7
-        self.assertEqual(len(new_reaction.genes), 1)
-        self.assertEqual(new_reaction.gene_reaction_rule, original_PGI_gpr)
-        # what about combining 2 gpr's
-        model.reactions.ACKr += model.reactions.ACONTa
-        self.assertEqual(model.reactions.ACKr.gene_reaction_rule,
-                         "(b2296 or b3115 or b1849) and (b0118 or b1276)")
-        self.assertEqual(len(model.reactions.ACKr.genes), 5)
-
-    def test_add(self):
-        # not in place addition should work on a copy
-        model = self.model
-        new = model.reactions.PGI + model.reactions.EX_h2o_e
-        self.assertIsNot(new._model, model)
-        self.assertEqual(len(new.metabolites), 3)
-        # the copy should refer to different metabolites and genes
-        # This currently fails because add_metabolites does not copy.
-        # Should that be changed?
-        # for met in new.metabolites:
-        #    self.assertIsNot(met, model.metabolites.get_by_id(met.id))
-        #    self.assertIsNot(met.model, model)
-        for gene in new.genes:
-            self.assertIsNot(gene, model.genes.get_by_id(gene.id))
-            self.assertIsNot(gene.model, model)
-
-    def test_mul(self):
-        new = self.model.reactions.PGI * 2
-        self.assertEqual(set(new.metabolites.values()), {-2, 2})
-
-    def test_sub(self):
-        model = self.model
-        new = model.reactions.PGI - model.reactions.EX_h2o_e
-        self.assertIsNot(new._model, model)
-        self.assertEqual(len(new.metabolites), 3)
-
-
-class TestCobraMetabolites(CobraTestCase):
-    def test_metabolite_formula(self):
-        met = Metabolite("water")
-        met.formula = "H2O"
-        self.assertEqual(met.elements, {"H": 2, "O": 1})
-        self.assertEqual(met.formula_weight, 18.01528)
-
-    def test_foruma_element_setting(self):
-        model = self.model
-        met = model.metabolites[1]
-        orig_formula = str(met.formula)
-        orig_elements = dict(met.elements)
-
-        met.formula = ''
-        self.assertEqual(met.elements, {})
-        met.elements = orig_elements
-        self.assertEqual(met.formula, orig_formula)
-
-
-class TestCobraModel(CobraTestCase):
-    """test core cobra functions"""
-
-    def test_add_reaction(self):
-        old_reaction_count = len(self.model.reactions)
-        old_metabolite_count = len(self.model.metabolites)
-        dummy_metabolite_1 = Metabolite("test_foo_1")
-        dummy_metabolite_2 = Metabolite("test_foo_2")
-        actual_metabolite = self.model.metabolites[0]
-        copy_metabolite = self.model.metabolites[1].copy()
-        dummy_reaction = Reaction("test_foo_reaction")
-        dummy_reaction.add_metabolites({dummy_metabolite_1: -1,
-                                        dummy_metabolite_2: 1,
-                                        copy_metabolite: -2,
-                                        actual_metabolite: 1})
-        self.model.add_reaction(dummy_reaction)
-        self.assertEqual(self.model.reactions.get_by_id(dummy_reaction.id),
-                         dummy_reaction)
-        for x in [dummy_metabolite_1, dummy_metabolite_2]:
-            self.assertEqual(self.model.metabolites.get_by_id(x.id), x)
-        # should have added 1 reaction and 2 metabolites
-        self.assertEqual(len(self.model.reactions), old_reaction_count + 1)
-        self.assertEqual(len(self.model.metabolites), old_metabolite_count + 2)
-        # tests on theadded reaction
-        reaction_in_model = self.model.reactions.get_by_id(dummy_reaction.id)
-        self.assertIs(type(reaction_in_model), Reaction)
-        self.assertIs(reaction_in_model, dummy_reaction)
-        self.assertEqual(len(reaction_in_model._metabolites), 4)
-        for i in reaction_in_model._metabolites:
-            self.assertEqual(type(i), Metabolite)
-        # tests on the added metabolites
-        met1_in_model = self.model.metabolites.get_by_id(dummy_metabolite_1.id)
-        self.assertIs(met1_in_model, dummy_metabolite_1)
-        copy_in_model = self.model.metabolites.get_by_id(copy_metabolite.id)
-        self.assertIsNot(copy_metabolite, copy_in_model)
-        self.assertIs(type(copy_in_model), Metabolite)
-        self.assertTrue(dummy_reaction in actual_metabolite._reaction)
-        # test adding a different metabolite with the same name as an
-        # existing one uses the metabolite in the model
-        r2 = Reaction("test_foo_reaction2")
-        self.model.add_reaction(r2)
-        r2.add_metabolites({Metabolite(self.model.metabolites[0].id): 1})
-        self.assertIs(self.model.metabolites[0], list(r2._metabolites)[0])
-
-    def test_add_reaction_from_other_model(self):
-        model = self.model
-        other = model.copy()
-        for i in other.reactions:
-            i.id += "_other"
-        other.repair()
-        model.add_reactions(other.reactions)
-        # what if the other reaction has an error in its GPR
-        m1 = create_test_model("textbook")
-        m2 = create_test_model("textbook")
-        m1.reactions.PGI.remove_from_model()
-        m2.genes.b4025._reaction.clear()
-        m1.add_reaction(m2.reactions.PGI)
-
-    def test_model_remove_reaction(self):
-        old_reaction_count = len(self.model.reactions)
-        self.model.remove_reactions(["PGI"])
-        self.assertEqual(len(self.model.reactions), old_reaction_count - 1)
-        with self.assertRaises(KeyError):
-            self.model.reactions.get_by_id("PGI")
-        self.model.remove_reactions(self.model.reactions[:1])
-        self.assertEqual(len(self.model.reactions), old_reaction_count - 2)
-        tmp_metabolite = Metabolite("testing")
-        self.model.reactions[0].add_metabolites({tmp_metabolite: 1})
-        self.assertIn(tmp_metabolite, self.model.metabolites)
-        self.model.remove_reactions(self.model.reactions[:1],
-                                    remove_orphans=True)
-        self.assertNotIn(tmp_metabolite, self.model.metabolites)
-
-    def test_reaction_remove(self):
-        model = self.model
-        old_reaction_count = len(model.reactions)
-        tmp_metabolite = Metabolite("testing")
-        # Delete without removing orphan
-        model.reactions[0].add_metabolites({tmp_metabolite: 1})
-        self.assertEqual(len(tmp_metabolite.reactions), 1)
-        # esnsure the stoichiometry is still the same using different objects
-        removed_reaction = model.reactions[0]
-        original_stoich = {i.id: value for i, value
-                           in removed_reaction._metabolites.items()}
-        model.reactions[0].remove_from_model(remove_orphans=False)
-        self.assertEqual(len(original_stoich),
-                         len(removed_reaction._metabolites))
-        for met in removed_reaction._metabolites:
-            self.assertEqual(original_stoich[met.id],
-                             removed_reaction._metabolites[met])
-            self.assertIsNot(met, model.metabolites)
-        # make sure it's still in the model
-        self.assertIn(tmp_metabolite, model.metabolites)
-        self.assertEqual(len(tmp_metabolite.reactions), 0)
-        self.assertEqual(len(self.model.reactions), old_reaction_count - 1)
-
-        # Now try it with removing orphans
-        model.reactions[0].add_metabolites({tmp_metabolite: 1})
-        self.assertEqual(len(tmp_metabolite.reactions), 1)
-        model.reactions[0].remove_from_model(remove_orphans=True)
-        self.assertNotIn(tmp_metabolite, model.metabolites)
-        self.assertEqual(len(tmp_metabolite.reactions), 0)
-        self.assertEqual(len(self.model.reactions), old_reaction_count - 2)
-
-        # It shouldn't remove orphans if it's in 2 reactions however
-        model.reactions[0].add_metabolites({tmp_metabolite: 1})
-        model.reactions[1].add_metabolites({tmp_metabolite: 1})
-        self.assertEqual(len(tmp_metabolite.reactions), 2)
-        model.reactions[0].remove_from_model(remove_orphans=False)
-        self.assertIn(tmp_metabolite, model.metabolites)
-        self.assertEqual(len(tmp_metabolite.reactions), 1)
-        self.assertEqual(len(self.model.reactions), old_reaction_count - 3)
-
-    def test_reaction_delete(self):
-        model = self.model
-        old_reaction_count = len(model.reactions)
-        tmp_metabolite = Metabolite("testing")
-        # Delete without removing orphan
-        model.reactions[0].add_metabolites({tmp_metabolite: 1})
-        self.assertEqual(len(tmp_metabolite.reactions), 1)
-        model.reactions[0].delete(remove_orphans=False)
-        # make sure it's still in the model
-        self.assertIn(tmp_metabolite, model.metabolites)
-        self.assertEqual(len(tmp_metabolite.reactions), 0)
-        self.assertEqual(len(self.model.reactions), old_reaction_count - 1)
-
-        # Now try it with removing orphans
-        model.reactions[0].add_metabolites({tmp_metabolite: 1})
-        self.assertEqual(len(tmp_metabolite.reactions), 1)
-        model.reactions[0].delete(remove_orphans=True)
-        self.assertNotIn(tmp_metabolite, model.metabolites)
-        self.assertEqual(len(tmp_metabolite.reactions), 0)
-        self.assertEqual(len(self.model.reactions), old_reaction_count - 2)
-
-        # It shouldn't remove orphans if it's in 2 reactions however
-        model.reactions[0].add_metabolites({tmp_metabolite: 1})
-        model.reactions[1].add_metabolites({tmp_metabolite: 1})
-        self.assertEqual(len(tmp_metabolite.reactions), 2)
-        model.reactions[0].delete(remove_orphans=False)
-        self.assertIn(tmp_metabolite, model.metabolites)
-        self.assertEqual(len(tmp_metabolite.reactions), 1)
-        self.assertEqual(len(self.model.reactions), old_reaction_count - 3)
-
-    def test_remove_gene(self):
-        target_gene = self.model.genes[0]
-        gene_reactions = list(target_gene.reactions)
-        with warnings.catch_warnings():
-            warnings.simplefilter("ignore")
-            target_gene.remove_from_model()
-        self.assertEqual(target_gene.model, None)
-        # make sure the reaction was removed from the model
-        self.assertNotIn(target_gene, self.model.genes)
-        # ensure the old reactions no longer have a record of the gene
-        for reaction in gene_reactions:
-            self.assertNotIn(target_gene, reaction.genes)
-
-    def test_copy(self):
-        """modifying copy should not modify the original"""
-        # test that deleting reactions in the copy does not change the
-        # number of reactions in the original model
-        model_copy = self.model.copy()
-        old_reaction_count = len(self.model.reactions)
-        self.assertEqual(len(self.model.reactions), len(model_copy.reactions))
-        self.assertEqual(len(self.model.metabolites),
-                         len(model_copy.metabolites))
-        model_copy.remove_reactions(model_copy.reactions[0:5])
-        self.assertEqual(old_reaction_count, len(self.model.reactions))
-        self.assertNotEqual(len(self.model.reactions),
-                            len(model_copy.reactions))
-
-    def test_deepcopy(self):
-        """Reference structures are maintained when deepcopying"""
-        model_copy = deepcopy(self.model)
-        for gene, gene_copy in zip(self.model.genes, model_copy.genes):
-            self.assertEqual(gene.id, gene_copy.id)
-            reactions = sorted(i.id for i in gene.reactions)
-            reactions_copy = sorted(i.id for i in gene_copy.reactions)
-            self.assertEqual(reactions, reactions_copy)
-        for reaction, reaction_copy in zip(self.model.reactions,
-                                           model_copy.reactions):
-            self.assertEqual(reaction.id, reaction_copy.id)
-            metabolites = sorted(i.id for i in reaction._metabolites)
-            metabolites_copy = sorted(i.id for i in reaction_copy._metabolites)
-            self.assertEqual(metabolites, metabolites_copy)
-
-    def test_add_reaction_orphans(self):
-        """test reaction addition
-
-        Need to verify that no orphan genes or metabolites are
-        contained in reactions after adding them to the model.
-        """
-        _model = self.model_class('test')
-        _model.add_reactions((x.copy() for x in self.model.reactions))
-        _genes = []
-        _metabolites = []
-        for x in _model.reactions:
-            _genes.extend(x.genes)
-            _metabolites.extend(x._metabolites)
-
-        orphan_genes = [x for x in _genes if x.model is not _model]
-        orphan_metabolites = [x for x in _metabolites if x.model is not _model]
-        self.assertEqual(len(orphan_genes), 0,
-                         msg='It looks like there are dangling genes when '
-                         'running Model.add_reactions')
-        self.assertEqual(len(orphan_metabolites), 0,
-                         msg='It looks like there are dangling metabolites '
-                         'when running Model.add_reactions')
-
-    def test_change_objective(self):
-        biomass = self.model.reactions.get_by_id("Biomass_Ecoli_core")
-        atpm = self.model.reactions.get_by_id("ATPM")
-        self.model.objective = atpm.id
-        self.assertEqual(atpm.objective_coefficient, 1.)
-        self.assertEqual(biomass.objective_coefficient, 0.)
-        self.assertEqual(self.model.objective, {atpm: 1.})
-        # change it back using object itself
-        self.model.objective = biomass
-        self.assertEqual(atpm.objective_coefficient, 0.)
-        self.assertEqual(biomass.objective_coefficient, 1.)
-        # set both to 1 with a list
-        self.model.objective = [atpm, biomass]
-        self.assertEqual(atpm.objective_coefficient, 1.)
-        self.assertEqual(biomass.objective_coefficient, 1.)
-        # set both using a dict
-        self.model.objective = {atpm: 0.2, biomass: 0.3}
-        self.assertEqual(atpm.objective_coefficient, 0.2)
-        self.assertEqual(biomass.objective_coefficient, 0.3)
-        # test setting by index
-        self.model.objective = self.model.reactions.index(atpm)
-        self.assertEqual(self.model.objective, {atpm: 1.})
-        # test by setting list of indexes
-        self.model.objective = map(self.model.reactions.index, [atpm, biomass])
-        self.assertEqual(self.model.objective, {atpm: 1., biomass: 1.})
-
-
- at skipIf(scipy is None, "scipy required for ArrayBasedModel")
-class TestCobraArrayModel(TestCobraModel):
-    def setUp(self):
-        model = create_test_model("textbook").to_array_based_model()
-        self.model_class = model.__class__
-        self.model = model
-
-    def test_array_based_model(self):
-        m = len(self.model.metabolites)
-        n = len(self.model.reactions)
-        assertEqual = self.assertEqual  # alias
-        for matrix_type in ["scipy.dok_matrix", "scipy.lil_matrix"]:
-            model = create_test_model("textbook").\
-                to_array_based_model(matrix_type=matrix_type)
-            assertEqual(model.S[7, 0], -1)
-            assertEqual(model.S[43, 0], 0)
-            model.S[43, 0] = 1
-            assertEqual(model.S[43, 0], 1)
-            assertEqual(
-                model.reactions[0]._metabolites[model.metabolites[43]], 1)
-            model.S[43, 0] = 0
-            assertEqual(model.lower_bounds[0], model.reactions[0].lower_bound)
-            assertEqual(model.lower_bounds[5], model.reactions[5].lower_bound)
-            assertEqual(model.upper_bounds[0], model.reactions[0].upper_bound)
-            assertEqual(model.upper_bounds[5], model.reactions[5].upper_bound)
-            model.lower_bounds[6] = 2
-            self.assertEqual(model.lower_bounds[6], 2)
-            self.assertEqual(model.reactions[6].lower_bound, 2)
-            # this should fail because it is the wrong size
-            with self.assertRaises(Exception):
-                model.upper_bounds = [0, 1]
-            model.upper_bounds = [0] * len(model.reactions)
-            self.assertEqual(max(model.upper_bounds), 0)
-
-            # test something for all the attributes
-            model.lower_bounds[2] = -1
-            assertEqual(model.reactions[2].lower_bound, -1)
-            assertEqual(model.lower_bounds[2], -1)
-            model.objective_coefficients[2] = 1
-            assertEqual(model.reactions[2].objective_coefficient, 1)
-            assertEqual(model.objective_coefficients[2], 1)
-            model.b[2] = 1
-            assertEqual(model.metabolites[2]._bound, 1)
-            assertEqual(model.b[2], 1)
-            model.constraint_sense[2] = "L"
-            assertEqual(model.metabolites[2]._constraint_sense, "L")
-            assertEqual(model.constraint_sense[2], "L")
-
-            # test resize matrix on reaction removal
-            m, n = model.S.shape
-            model.remove_reactions([model.reactions[2]], remove_orphans=False)
-            self.assertEqual(len(model.metabolites), model.S.shape[0])
-            self.assertEqual(len(model.reactions), model.S.shape[1])
-            self.assertEqual(model.S.shape, (m, n - 1))
-
-    def test_array_based_model_add(self):
-        m = len(self.model.metabolites)
-        n = len(self.model.reactions)
-        for matrix_type in ["scipy.dok_matrix", "scipy.lil_matrix"]:
-            model = create_test_model("textbook").\
-                to_array_based_model(matrix_type=matrix_type)
-            test_reaction = Reaction("test")
-            test_reaction.add_metabolites({model.metabolites[0]: 4})
-            test_reaction.lower_bound = -3.14
-            model.add_reaction(test_reaction)
-            self.assertEqual(len(model.reactions), n + 1)
-            self.assertEqual(model.S.shape, (m, n + 1))
-            self.assertEqual(len(model.lower_bounds), n + 1)
-            self.assertEqual(len(model.upper_bounds), n + 1)
-            self.assertEqual(model.S[0, n], 4)
-            self.assertEqual(model.S[7, 0], -1)
-            self.assertEqual(model.lower_bounds[n], -3.14)
-
-    def test_array_based_select(self):
-        model = self.model
-        atpm_select = model.reactions[model.lower_bounds > 0]
-        self.assertEqual(len(atpm_select), 1)
-        self.assertEqual(atpm_select[0].id, "ATPM")
-        self.assertEqual(len(model.reactions[model.lower_bounds <= 0]),
-                         len(model.reactions) - 1)
-        # mismatched dimensions should give an error
-        with self.assertRaises(TypeError):
-            model.reactions[[True, False]]
-
-    def test_array_based_bounds_setting(self):
-        model = self.model
-        bounds = [0.0] * len(model.reactions)
-        model.lower_bounds = bounds
-        self.assertEqual(type(model.reactions[0].lower_bound), float)
-        self.assertAlmostEqual(model.reactions[0].lower_bound, 0.0)
-        model.upper_bounds[1] = 1234.0
-        self.assertAlmostEqual(model.reactions[1].upper_bound, 1234.0)
-        model.upper_bounds[9:11] = [100.0, 200.0]
-        self.assertAlmostEqual(model.reactions[9].upper_bound, 100.0)
-        self.assertAlmostEqual(model.reactions[10].upper_bound, 200.0)
-        model.upper_bounds[9:11] = 123.0
-        self.assertAlmostEqual(model.reactions[9].upper_bound, 123.0)
-        self.assertAlmostEqual(model.reactions[10].upper_bound, 123.0)
-
-
-# make a test suite to run all of the tests
-loader = TestLoader()
-suite = loader.loadTestsFromModule(sys.modules[__name__])
-
-
-def test_all():
-    TextTestRunner(verbosity=2).run(suite)
-
-if __name__ == "__main__":
-    test_all()
diff --git a/config.sh b/config.sh
index b8ba8a7..b855e0a 100644
--- a/config.sh
+++ b/config.sh
@@ -7,26 +7,25 @@ function pre_build {
     if [ -n "$IS_OSX" ]; then
         export CC=clang
         export CXX=clang++
-        export CFLAGS="-fPIC -O3 -arch i386 -arch x86_64 -g -DNDEBUG -mmacosx-version-min=10.6"
-    else
-        yum install -y libxslt libxml2 libxml2-devel libxslt-devel
-    fi
-    curl -O http://ftp.gnu.org/gnu/glpk/glpk-4.60.tar.gz
-    tar xzf glpk-4.60.tar.gz
-    (cd glpk-4.60 \
-            && ./configure \
-            && make \
-            && make install)
-    pip install cython
-    cython -a cobra/solvers/cglpk.pyx
-    export PATH="$PATH:/usr/local/bin"
+		export CFLAGS="-fPIC -O3 -arch i386 -arch x86_64 -g -DNDEBUG -mmacosx-version-min=10.6"
+	else
+		yum install -y libxslt libxml2 libxml2-devel libxslt-devel
+	fi
+	curl -O http://ftp.gnu.org/gnu/glpk/glpk-4.60.tar.gz
+	tar xzf glpk-4.60.tar.gz
+	(cd glpk-4.60 \
+			&& ./configure \
+			&& make \
+			&& make install)
+	pip install cython
+	cython -a cobra/solvers/cglpk.pyx
+	export PATH="$PATH:/usr/local/bin"
 }
 
 function build_wheel {
     # Set default building method to pip
     build_bdist_wheel $@
-	# avoid this for now, (we get broken linux wheels due to https://github.com/pypa/manylinux/issues/80, but testing works)
-	#(cd glpk-4.60 && make uninstall)
+	(cd glpk-4.60 && make uninstall)
 }
 
 function run_tests_in_repo {
@@ -46,9 +45,7 @@ function run_tests_in_repo {
 	fi
 	mkdir -p $HOME/.config/matplotlib
 	echo 'backend: Agg' >> $HOME/.config/matplotlib/matplotlibrc
-	echo -e "import cobra.test; import sys; sys.exit(cobra.test.test_all())" > run-tests.py
-	(coverage run --source=cobra --rcfile ../.coveragerc run-tests.py &&
-			coverage xml &&
+	(pytest --pyargs -v -rsx --cov=cobra --cov-report=xml --cov-config=../.coveragerc --benchmark-skip cobra &&
 			mv coverage.xml ..)
 }
 
diff --git a/develop-requirements.txt b/develop-requirements.txt
index f72cca2..2380179 100644
--- a/develop-requirements.txt
+++ b/develop-requirements.txt
@@ -13,3 +13,4 @@ pandas>=0.17.0
 tabulate
 tox
 pep8
+pytest
diff --git a/documentation_builder/conf.py b/documentation_builder/conf.py
index ae009a1..bc80bdf 100644
--- a/documentation_builder/conf.py
+++ b/documentation_builder/conf.py
@@ -37,7 +37,7 @@ class Mock(object):
 
 MOCK_MODULES = ['numpy', 'scipy', 'scipy.sparse', 'scipy.io', 'scipy.stats',
                 'glpk', 'gurobipy', 'gurobipy.GRB', 'cplex', 'pp', 'libsbml',
-                'cplex.exceptions', 'pandas']
+                'cplex.exceptions', 'pandas', 'tabulate']
 for mod_name in MOCK_MODULES:
     sys.modules[mod_name] = Mock()
 
diff --git a/scripts/deploy-test.sh b/scripts/deploy-test.sh
deleted file mode 100755
index 6db962f..0000000
--- a/scripts/deploy-test.sh
+++ /dev/null
@@ -1,15 +0,0 @@
-#!/bin/bash
-echo -e " ... running twine to deploy ... "
-echo -e "
-[distutils]
-index-servers=
-    pypirepository
-
-[pypirepository]
-repository: https://testpypi.python.org/pypi
-username: henred
-password: pippi_langstrump
-" > ~/.pypirc
-
-pip install twine
-twine upload --skip-existing ${TRAVIS_BUILD_DIR}/wheelhouse/* -r pypirepository
diff --git a/scripts/deploy.sh b/scripts/deploy.sh
index fe46af5..fcebd2f 100755
--- a/scripts/deploy.sh
+++ b/scripts/deploy.sh
@@ -1,6 +1,7 @@
 #!/bin/bash
-echo -e " ... running twine to deploy ... "
-pip install twine
 
-# only upload mac wheels for now, linux wheels are broken due to https://github.com/pypa/manylinux/issues/80
-twine upload --skip-existing --username "${PYPI_USERNAME}" --password "${PYPI_PASSWORD}" ${TRAVIS_BUILD_DIR}/wheelhouse/*macosx*.whl
+
+echo -e " starting deploy for branch ${TRAVIS_BRANCH} .."
+
+pip install twine
+twine upload --skip-existing --username "${PYPI_USERNAME}" --password "${PYPI_PASSWORD}" ${TRAVIS_BUILD_DIR}/wheelhouse/*
diff --git a/tox.ini b/tox.ini
index 6c4f934..68ef7fa 100644
--- a/tox.ini
+++ b/tox.ini
@@ -9,4 +9,4 @@ commands=pep8 cobra --exclude=oven,solvers,sbml.py
 [testenv]
 commands =
     pip install -U pip
-    python setup.py test
+    pytest --benchmark-skip

-- 
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/python-cobra.git



More information about the debian-med-commit mailing list