Fix issue #1300 (#1312)

* Fix bug

* Added test

* Small fix from testing

* Made test a bit easier to understand

* Black reformatting

* Revert "Black reformatting"

This reverts commit 55116a7.

* Revert "Revert "Black reformatting""

This reverts commit 4c6a16e.

* Added next release notes

release-notes/next-release.md

@@ -3,6 +3,7 @@
 ## New features
 
 ## Fixes
+[#1300](https://github.com/opencobra/cobrapy/pull/1312)- Fixed issue where reaction bounds were being overwritten by global model min/max values when writing sbml model to file
 
 ## Other
 

src/cobra/core/formula.py

@@ -63,7 +63,7 @@ def parse_composition(self) -> None:
             return
         composition = {}
         parsed = element_re.findall(tmp_formula)
-        for (element, count) in parsed:
+        for element, count in parsed:
             if count == "":
                 count = 1
             else:

src/cobra/core/metabolite.py

@@ -133,7 +133,7 @@ def elements(self) -> Optional[Dict[str, Union[int, float]]]:
             return None
         composition = {}
         parsed = element_re.findall(tmp_formula)
-        for (element, count) in parsed:
+        for element, count in parsed:
             if count == "":
                 count = 1
             else:

src/cobra/core/model.py

@@ -32,7 +32,6 @@
 
 
 if TYPE_CHECKING:
-
     import pandas as pd
     from optlang.container import Container
 
@@ -781,7 +780,6 @@ def remove_reactions(
         context = get_context(self)
 
         for reaction in reactions:
-
             # Make sure the reaction is in the model
             try:
                 reaction = self.reactions[self.reactions.index(reaction)]
@@ -793,7 +791,6 @@ def remove_reactions(
                 reverse = reaction.reverse_variable
 
                 if context:
-
                     obj_coef = reaction.objective_coefficient
 
                     if obj_coef != 0:

src/cobra/core/reaction.py

@@ -1228,7 +1228,6 @@ def add_metabolites(
         _id_to_metabolites = dict([(x.id, x) for x in self._metabolites])
 
         for metabolite, coefficient in metabolites_to_add.items():
-
             # Make sure metabolites being added belong to the same model, or
             # else copy them.
             if isinstance(metabolite, Metabolite):

src/cobra/core/solution.py

@@ -177,21 +177,21 @@ def get_solution(
     if model.solver.is_integer:
         reduced.fill(np.nan)
         shadow.fill(np.nan)
-        for (i, rxn) in enumerate(reactions):
+        for i, rxn in enumerate(reactions):
             rxn_index.append(rxn.id)
             fluxes[i] = var_primals[rxn.id] - var_primals[rxn.reverse_id]
         met_index = [met.id for met in metabolites]
     else:
         var_duals = model.solver.reduced_costs
-        for (i, rxn) in enumerate(reactions):
+        for i, rxn in enumerate(reactions):
             forward = rxn.id
             reverse = rxn.reverse_id
             rxn_index.append(forward)
             fluxes[i] = var_primals[forward] - var_primals[reverse]
             reduced[i] = var_duals[forward] - var_duals[reverse]
         met_index = []
         constr_duals = model.solver.shadow_prices
-        for (i, met) in enumerate(metabolites):
+        for i, met in enumerate(metabolites):
             met_index.append(met.id)
             shadow[i] = constr_duals[met.id]
     return Solution(

src/cobra/io/sbml.py

@@ -1038,7 +1038,6 @@ def process_association(ass: "libsbml.FbcAssociation") -> Union[BoolOp, Name]:
             model.getListOfReactions(),
             model_groups.getListOfGroups(),
         ]:
-
             sbase: "libsbml.SBase"
             for sbase in obj_list:
                 if sbase.isSetId():
@@ -1277,13 +1276,8 @@ def _model_to_sbml(
             unit.setScale(u.scale)
             unit.setMultiplier(u.multiplier)
 
-    # minimum and maximum value from model
-    if len(cobra_model.reactions) > 0:
-        min_value = min(cobra_model.reactions.list_attr("lower_bound"))
-        max_value = max(cobra_model.reactions.list_attr("upper_bound"))
-    else:
-        min_value = config.lower_bound
-        max_value = config.upper_bound
+    min_value = config.lower_bound
+    max_value = config.upper_bound
 
     _create_parameter(
         model, pid=LOWER_BOUND_ID, value=min_value, sbo=SBO_DEFAULT_FLUX_BOUND
@@ -1875,7 +1869,6 @@ def _sbase_annotations(sbase: libsbml.SBase, annotation: dict) -> None:
 
     # rdf_items = []
     for provider, data in annotation_data.items():
-
         # set SBOTerm
         if provider in ["SBO", "sbo"]:
             if provider == "SBO":

src/cobra/sampling/hr_sampler.py

@@ -323,9 +323,7 @@ def generate_fva_warmup(self) -> None:
 
                 primals = self.model.solver.primal_values
                 sol = [primals[v.name] for v in self.model.variables]
-                self.warmup[
-                    self.n_warmup,
-                ] = sol
+                self.warmup[self.n_warmup,] = sol
                 self.n_warmup += 1
 
                 # Reset objective
@@ -425,33 +423,13 @@ def _is_redundant(self, matrix: np.matrix, cutoff: Optional[float] = None) -> bo
     def _bounds_dist(self, p: np.ndarray) -> np.ndarray:
         """Get the lower and upper bound distances. Negative is bad."""
         prob = self.problem
-        lb_dist = (
-            p
-            - prob.variable_bounds[
-                0,
-            ]
-        ).min()
-        ub_dist = (
-            prob.variable_bounds[
-                1,
-            ]
-            - p
-        ).min()
+        lb_dist = (p - prob.variable_bounds[0,]).min()
+        ub_dist = (prob.variable_bounds[1,] - p).min()
 
         if prob.bounds.shape[0] > 0:
             const = prob.inequalities.dot(p)
-            const_lb_dist = (
-                const
-                - prob.bounds[
-                    0,
-                ]
-            ).min()
-            const_ub_dist = (
-                prob.bounds[
-                    1,
-                ]
-                - const
-            ).min()
+            const_lb_dist = (const - prob.bounds[0,]).min()
+            const_ub_dist = (prob.bounds[1,] - const).min()
             lb_dist = min(lb_dist, const_lb_dist)
             ub_dist = min(ub_dist, const_ub_dist)
 
@@ -564,38 +542,18 @@ def validate(self, samples: np.matrix) -> np.ndarray:
             )
 
         feasibility = np.abs(S.dot(samples.T).T - b).max(axis=1)
-        lb_error = (
-            samples
-            - bounds[
-                0,
-            ]
-        ).min(axis=1)
-        ub_error = (
-            bounds[
-                1,
-            ]
-            - samples
-        ).min(axis=1)
+        lb_error = (samples - bounds[0,]).min(axis=1)
+        ub_error = (bounds[1,] - samples).min(axis=1)
 
         if samples.shape[1] == len(self.model.variables) and prob.inequalities.shape[0]:
             consts = prob.inequalities.dot(samples.T)
             lb_error = np.minimum(
                 lb_error,
-                (
-                    consts
-                    - prob.bounds[
-                        0,
-                    ]
-                ).min(axis=1),
+                (consts - prob.bounds[0,]).min(axis=1),
             )
             ub_error = np.minimum(
                 ub_error,
-                (
-                    prob.bounds[
-                        1,
-                    ]
-                    - consts
-                ).min(axis=1),
+                (prob.bounds[1,] - consts).min(axis=1),
             )
 
         valid = (

tests/test_flux_analysis/test_deletion.py

@@ -428,7 +428,7 @@ def test_double_reaction_deletion(model: Model) -> None:
     solution_one_process = double_reaction_deletion(
         model, reaction_list1=reactions, processes=1
     )
-    for (rxn_a, sub) in growth_dict.items():
+    for rxn_a, sub in growth_dict.items():
         for rxn_b, growth in sub.items():
             sol = solution.knockout[{rxn_a, rxn_b}]
             sol_one = solution_one_process.knockout[{rxn_a, rxn_b}]

tests/test_io/test_sbml.py

@@ -511,6 +511,45 @@ def test_boundary_conditions(data_directory: Path) -> None:
     assert sol1.objective_value == sol2.objective_value
 
 
+def test_bounds_on_write(data_directory: Path, tmp_path: Path) -> None:
+    """Test infinity bound example.
+
+    Parameters
+    ----------
+    data_directory: Path
+        Directory where the data is.
+    """
+    sbml_path1 = data_directory / "fbc_ex1.xml"
+    model1 = read_sbml_model(sbml_path1)
+
+    r_x = model1.reactions.get_by_id("EX_X")
+    r_y = model1.reactions.get_by_id("EX_Ac")
+
+    r_x.bounds = (config.lower_bound - 1000, config.upper_bound + 1000)
+    assert r_x.lower_bound == config.lower_bound - 1000
+    assert r_x.upper_bound == config.upper_bound + 1000
+
+    # Global min/max bounds for other reactions should not change before & after write!
+    r_y.bounds = (config.lower_bound, config.upper_bound)
+    assert r_y.lower_bound == config.lower_bound
+    assert r_y.upper_bound == config.upper_bound
+
+    sbml_path = tmp_path / "test.xml"
+    with open(sbml_path, "w") as f_out:
+        write_sbml_model(model1, f_out)
+
+    with open(sbml_path, "r") as f_in:
+        model2 = read_sbml_model(f_in)
+
+    r2_x = model2.reactions.get_by_id("EX_X")
+    r2_y = model2.reactions.get_by_id("EX_Ac")
+
+    assert r2_x.lower_bound == r_x.lower_bound
+    assert r2_x.upper_bound == r_x.upper_bound
+    assert r2_y.lower_bound == r_y.lower_bound  # before fix #1300, this would fail
+    assert r2_y.upper_bound == r_y.upper_bound  # before fix #1300, this would fail
+
+
 def test_gprs(large_model: Model, tmp_path: Path) -> None:
     """Test that GPRs are written and read correctly.
 

tests/test_util/test_context.py

@@ -39,6 +39,7 @@ def test_get_context(model: "Model") -> None:
 
 def test_resettable() -> None:
     """Test if resettable decorator is functional."""
+
     # decorate a dummy function
     @resettable
     def change_my_name(old_name, new_name):