fix(chemistry): fix all negative elimination coefs

Catch and fix the case of all negative elimination coefficients during
elimination.  Improve documentation and comments in the equilibrium
eliminate method.  Add, parametrize, and clarify tests for the
equilibrium eliminate method.

Github-closes: closes bjodah#209
Signed-off-by: Jeremy A Gray <jeremy.a.gray@gmail.com>

chempy/chemistry.py

@@ -1237,7 +1237,11 @@ def __sub__(self, other):
 
     @staticmethod
     def eliminate(rxns, wrt):
-        """Linear combination coefficients for elimination of a substance
+        """Eliminate a substance from a pair of reactions.
+
+        Find the linear combination coefficients for elimination of a
+        substance in a pair of reactions.  Reactions should be
+        balanced before elimination.
 
         Parameters
         ----------
@@ -1256,15 +1260,35 @@ def eliminate(rxns, wrt):
         """
         import sympy
 
+        # Reactants have negative coefficients; products positive.
         viol = [r.net_stoich([wrt])[0] for r in rxns]
         factors = defaultdict(int)
+
+        # Get the prime factorization for each coefficient and keep
+        # larger count of factors.
         for v in viol:
-            for k, v in sympy.factorint(v).items():
-                factors[k] = max(factors[k], v)
+            for key, val in sympy.factorint(v).items():
+                factors[key] = max(factors[key], val)
+
+        # Calculate the LCM, which is the product of the species
+        # coefficient and the elimination coefficient for the
+        # reaction.
         rcd = reduce(mul, (k ** v for k, v in factors.items()))
+
+        # Negate one coefficient to eliminate (subtract).  The
+        # coefficient check below will catch any all negative sets of
+        # coefficients.
         viol[0] *= -1
 
-        return [rcd // v for v in viol]
+        # Divide the LCM by each species coeffient to calculate the
+        # elimination coefficient for each reaction.
+        coefs = [rcd // v for v in viol]
+
+        # If all the coefficients are non-positive, negate them to
+        # avoid reversing the reaction unnecessarily.
+        coefs = [-v for v in coefs] if all([v <= 0 for v in coefs]) else coefs
+
+        return coefs
 
     def cancel(self, rxn):
         """Multiplier of how many times rxn can be added/subtracted.

chempy/tests/test_chemistry.py

@@ -303,6 +303,28 @@ def test_Equilibrium__eliminate():
     assert (e2 * -3).reac == {"E": 33}
 
 
+@pytest.mark.parametrize(
+    "r1, p1, r2, p2, wrt, c",
+    [
+        ({"A": 1, "B": 2}, {"C": 3}, {"D": 5, "B": 7}, {"E": 11}, "B", [-7, 2]),
+        (
+            {"Cd+2": 4, "H2O": 4},
+            {"Cd4(OH)4+4": 1, "H+": 4},
+            {"Cd(OH)2(s)": 1},
+            {"Cd+2": 1, "OH-": 2},
+            "Cd+2",
+            [1, 4],
+        ),
+    ],
+)
+@requires("sympy")
+def test_equilibrium_eliminate(r1, p1, r2, p2, wrt, c):
+    e1 = Equilibrium(r1, p1)
+    e2 = Equilibrium(r2, p2)
+    coeff = Equilibrium.eliminate([e1, e2], wrt)
+    assert c == coeff
+
+
 @pytest.mark.parametrize(
     "r1, p1, r2, p2, wrt, c",
     [
@@ -320,7 +342,7 @@ def test_Equilibrium__eliminate():
             {"Br-", "H2O"},
             {"BrO3-", "H+", "e-"},
             "e-",
-            [-2, -1],
+            [2, 1],
         ),
     ],
 )
@@ -334,7 +356,7 @@ def test_equilibrium_balance_eliminate_coefficients(r1, p1, r2, p2, wrt, c):
 
     coeff = Equilibrium.eliminate([one, two], wrt)
 
-    assert coeff == c
+    assert c == coeff
 
 
 @requires(parsing_library, units_library)