Add custom solution class, several fixes

- adds a custom solution class
- more tests
- optimize will now use pfba by default if fluxes are requested
- slim option to only get growth rates

micom/solution.py

@@ -0,0 +1,89 @@
+"""A community solution object."""
+
+import numpy as np
+import pandas as pd
+from optlang.interface import OPTIMAL
+from cobra.core import Solution, get_solution
+
+
+def _group_species(values, ids, species, what="reaction"):
+    """Format a list of values by id and species."""
+    df = pd.DataFrame({values.name: values, what: ids, "species": species})
+    df = df.pivot(index="species", columns=what, values=values.name)
+    df.name = values.name
+    return df
+
+
+class CommunitySolution(Solution):
+    """An FBA solution for an entire community.
+
+    Attributes
+    ----------
+    objective_value : float
+        The (optimal) value for the objective function.
+    growth_rates : pandas.Series
+        The growth rates for each of the species.
+    community_growth : float
+        The overall growth rate for the community normalized to 1 gDW.
+    status : str
+        The solver status related to the solution.
+    reactions : numpy.array
+        A list of `cobra.Reaction` objects for which the solution is
+        retrieved.
+    fluxes : pandas.DataFrame
+        Contains the reaction fluxes (primal values of variables) stratified
+        by species. Columns denote individual fluxes and rows denote species.
+        Fluxes will be NA if the reaction does not exist in the organism.
+    reduced_costs : pandas.Series
+        Contains reaction reduced costs (dual values of variables) stratified
+        by species. Columns denote individual fluxes and rows denote species.
+        Reduced costs will be NA if the reaction does not exist in the
+        organism.
+    metabolites : numpy.array
+        A list of `cobra.Metabolite` objects for which the solution is
+        retrieved.
+    shadow_prices : pandas.Series
+        Contains metabolite shadow prices (dual values of constraints)
+        stratified by species. Columns denote individual metabolites and rows
+        denote species. Shadow prices will be NA if the metabolite does not
+        exist in the organism.
+
+    """
+
+    def __init__(self, community, slim=False,
+                 reactions=None, metabolites=None):
+        """Get the solution from a community model."""
+        if reactions is None:
+            reactions = community.reactions
+        if metabolites is None:
+            metabolites = community.metabolites
+        if not slim:
+            rids = np.array([(r.global_id, r.community_id) for r in reactions])
+            mids = np.array([(m.global_id, m.community_id)
+                             for m in metabolites])
+            sol = get_solution(community, reactions, metabolites)
+            super(CommunitySolution, self).__init__(
+                community.solver.objective.value, community.solver.status,
+                np.unique(rids[:, 0]),
+                _group_species(sol.fluxes, rids[:, 0], rids[:, 1]),
+                _group_species(sol.reduced_costs, rids[:, 0], rids[:, 1]),
+                np.unique(mids[:, 0]),
+                _group_species(sol.shadow_prices, mids[:, 0], mids[:, 1],
+                               what="metabolites"))
+        else:
+            super(CommunitySolution, self).__init__(
+                community.solver.objective.value, community.solver.status,
+                None, None, None, None, None)
+        gcs = pd.Series()
+        for sp in community.objectives:
+            gcs[sp] = community.constraints["objective_" + sp].primal
+        self.growth_rates = gcs
+        self.community_growth = sum(community.abundances * gcs)
+
+    def __repr__(self):
+        """Convert CommunitySolution instance to string representation."""
+        if self.status != OPTIMAL:
+            return "<CommunitySolution {0:s} at 0x{1:x}>".format(
+                self.status, id(self))
+        return "<CommunitySolution {0:.3f} at 0x{1:x}>".format(
+            self.community_growth, id(self))

tests/test_media.py

@@ -0,0 +1,31 @@
+"""Test growth media agorithms."""
+
+from fixtures import community
+import numpy as np
+import micom.media as media
+
+
+class TestMinimalMedia:
+
+    def test_linear(self, community):
+        medium = media.minimal_medium(community, 0.8, 0.1)
+        assert len(medium) <= 4
+        assert all(medium > 1e-9)
+
+    def test_mip(self, community):
+        medium = media.minimal_medium(community, 0.8, 0.1,
+                                      minimize_components=True)
+        assert len(medium) <= 4
+        assert all(medium > 1e-9)
+
+        # Anaerobic growth
+        medium = media.minimal_medium(community, 0.1, 0.1,
+                                      minimize_components=True)
+        assert len(medium) <= 3
+        assert all(medium > 1e-9)
+
+    def test_benchmark_linear(self, community, benchmark):
+        benchmark(media.minimal_medium, community, 0.8, 0.1)
+
+    def test_benchmark_mip(self, community, benchmark):
+        benchmark(media.minimal_medium, community, 0.8, 0.1, True)