Skip to content

Commit

Permalink
Merge pull request #58 from CalebBell/caleb/difficult_reaction_balancing
Browse files Browse the repository at this point in the history 
Caleb/difficult reaction balancing
  • Loading branch information
@CalebBell
CalebBell authored Nov 20, 2024
2 parents 81476a4 + 1ac5220 commit 6e6d3b4
Show file tree
Hide file tree
Showing 2 changed files with 150 additions and 58 deletions.
85 changes: 48 additions & 37 deletions chemicals/reaction.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -995,11 +995,42 @@ def stoichiometric_matrix(atomss, reactants):
def round_to_significant(x, significant_digits):
if x == 0:
return 0.0

magnitude = floor(log10(abs(x)))
scale = 10 ** (significant_digits - 1 - magnitude)
return round(x * scale) / scale

def check_reaction_balance(matrix, coeffs, atol=1e-13):
"""Check that coefficients satisfy the stoichiometric matrix equation within tolerance."""
result = [sum(coeff * row[i] for i, coeff in enumerate(coeffs))
for row in matrix]
return all(abs(x) <= atol for x in result)

def floats_to_ints(float_list, matrix, max_denominator=1000):
"""
Convert a list of floats to integers until we find a solution that balances.
All floats are one or larger. The chemical equation is assumed to be reasonable.
The SVD has already solved the problem, but there is a little numerical noise
we need to clean up.
Parameters:
- float_list: List of floats to convert
- matrix: Stoichiometric matrix to verify balance
- max_denominator: Maximum scaling factor to consider
Returns:
- A list of integers scaled from the original floats
"""
for D in range(1, max_denominator + 1):
# Calculate rounded integers
# in practice this works extremely well and rarely goes above 10
# it is extremely fast compared to a Fraction/Decimal approach
rounded = [int(round(D * x)) for x in float_list]
# Check if these coefficients actually balance the reaction
if check_reaction_balance(matrix, rounded):
return rounded
return float_list # If we still can't find a solution, return original floats


def balance_stoichiometry(matrix, rounding=9, allow_fractional=False):
r'''This function balances a chemical reaction.
Expand All @@ -1008,6 +1039,15 @@ def balance_stoichiometry(matrix, rounding=9, allow_fractional=False):
matrix : list[list[float]]
Chemical reaction matrix for further processing; rows contain element
counts of each compound, and the columns represent each chemical, [-]
rounding : int
Roughly the number of digits of rounding to apply to the answer. As matrix
routines are used, there is some noise; if this number is too high, the
coefficients may become very large numberes, which are still in a correct
ratio to each other, but are extremely ugly, [-]
allow_fractional : bool
Whether or not to allow the answers to be fractions, or to force them to
integers. Setting this to True speeds up the calculation, and allows
setting rounding arbitrarily high, [-]
Returns
-------
Expand All @@ -1033,6 +1073,7 @@ def balance_stoichiometry(matrix, rounding=9, allow_fractional=False):
This algorithm may suffer from floating point issues. If you believe there
is an error in the result, please report your reaction to the developers.
This function has a comprehensive test suite and extra test cases can be added to it.
References
----------
Expand All @@ -1050,47 +1091,17 @@ def balance_stoichiometry(matrix, rounding=9, allow_fractional=False):
74, no. 11 (November 1, 1997): 1369. https://doi.org/10.1021/ed074p1369.
'''
import scipy.linalg
done = scipy.linalg.null_space(matrix)
done = scipy.linalg.null_space(matrix, rcond=None)
if len(done[0]) > 1:
raise ValueError("No solution")
d = done[:, 0].tolist()

min_value_inv = 1.0/min(d)
min_value_inv = 1.0/min(d, key=abs)
d = [i*min_value_inv for i in d]

d = [round_to_significant(v, rounding) for v in d]
if not allow_fractional:
from fractions import Fraction
max_denominator = 10**rounding
fs = [Fraction(x).limit_denominator(max_denominator=max_denominator) for x in d]
all_denominators = {i.denominator for i in fs}
all_denominators.discard(Fraction(1))

for den in sorted(list(all_denominators), reverse=True):
fs = [num*den for num in fs]
if all(i.denominator == 1 for i in fs):
break

# May have gone too far
return [float(i) for i in fs]
# done = False
# for i in range(100):
# for c in d:
# ratio = c.as_integer_ratio()[1]
# if ratio != 1:
# d = [di*ratio for di in d]
# break
# done = True
# if done:
# break
#
# d_as_int = [int(i) for i in d]
# for i, j in zip(d, d_as_int):
# if i != j:
# raise ValueError("Could not find integer coefficients (%s, %s)" %(i, j))
# return d_as_int
else:
d = [round_to_significant(v, rounding) for v in d]
return d
d = floats_to_ints(d, matrix)
d = [float(v) for v in d]
return d

def stoichiometry_molar_to_mass(coefficients, MWs):
r'''This function translates molar stoichiometric
Expand Down
123 changes: 102 additions & 21 deletions tests/test_reaction.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -210,22 +210,97 @@ def check_reaction_balance(matrix, products_calc, atol=1e-13):
result = np.array(matrix) @ np.array(products_calc)
assert_close1d(result, [0.0]*len(result), atol=atol)

# def test_balance_stoichiometry_ill_conditioned():
# test_cases = [
# # C100000H200000N + O2 = CO2 + H2O + NO2
# [[{'C': 100000, 'H': 200000, 'N': 1}, {'O': 2}, {'C': 1, 'O': 2}, {'H': 2, 'O': 1}, {'N': 1, 'O': 2}],
# [True, True, False, False, False],
# [5.5322455442586566e+17, 8.298423638960436e+22, 5.532245544336663e+22, 5.532245544336663e+22, 5.5322455442586566e+17]],
# ]
# for atomss, statuses, products in test_cases:
# matrix = stoichiometric_matrix(atomss, statuses)
# products_calc = balance_stoichiometry(matrix)
# check_reaction_balance(matrix, products_calc)
# assert_close1d(products_calc, products)


def test_balance_stoichiometry():
test_cases = [
ill_conditioned_stoich_test_cases = [
(
[[{'C': 100000, 'H': 200000, 'N': 1}, {'O': 2}, {'C': 1, 'O': 2}, {'H': 2, 'O': 1}, {'N': 1, 'O': 2}],
[True, True, False, False, False],
[1, 150001, 100000, 100000, 1]],
"C100000H200000N + O2 = CO2 + H2O + NO2"
),
(
[[{'C': 100000000, 'H': 200000000, 'N': 1}, {'O': 2}, {'C': 1, 'O': 2}, {'H': 2, 'O': 1}, {'N': 1, 'O': 2}],
[True, True, False, False, False],
[1, 150000001, 100000000, 100000000, 1]],
"C100000000H200000000N + O2 = CO2 + H2O + NO2"
),
(
[[{'C': 50000, 'H': 100000, 'O': 25000}, {'O': 2}, {'C': 1, 'O': 2}, {'H': 2, 'O': 1}],
[True, True, False, False],
[1.0, 62500.0, 50000.0, 50000.0]],
"C50000H100000O25000 + O2 = CO2 + H2O"
),
(
[[{'Fe': 1000, 'C': 2000, 'N': 2000}, {'H': 2, 'O': 2}, {'Fe': 2, 'O': 3}, {'H': 1, 'C': 1, 'N': 1}, {'H': 2, 'O': 1}],
[True, True, False, False, False],
[1.0, 500.0, 500.0, 2000.0, -500.0]],
"Fe1000(CN)2000 + H2O2 = Fe2O3 + HCN + H2O"
),
(
[[{'Au': 1000, 'Cu': 1000}, {'H': 1, 'N': 1, 'O': 3}, {'Au': 1}, {'Cu': 1, 'N': 2, 'O': 6}, {'N': 1, 'O': 1}, {'H': 2, 'O': 1}],
[True, True, False, False, False, False],
[3, 8000, 3000, 3000, 2000, 4000]],
"Au1000Cu1000 + HNO3 = Au + Cu(NO3)2 + NO + H2O"
),
(
[[{'H': 100000, 'S': 50000, 'O': 200000}, {'S': 1, 'O': 3}, {'H': 2, 'O': 1}],
[True, False, False],
[1.0, 50000.0, 50000.0]],
"H100000(SO4)50000 = SO3 + H2O"
),
(
[[{'C': 50000, 'H': 100000, 'N': 10000, 'O': 20000, 'S': 1000}, {'O': 2}, {'C': 1, 'O': 2}, {'H': 2, 'O': 1}, {'N': 1, 'O': 2}, {'S': 1, 'O': 2}],
[True, True, False, False, False, False],
[1.0, 76000.0, 50000.0, 50000.0, 10000.0, 1000.0]],
"C50000H100000N10000O20000S1000 + O2 = CO2 + H2O + NO2 + SO2"
),
(
[[{'C': 12594, 'H': 25422, 'O': 5004}, {'O': 2}, {'C': 1, 'O': 2}, {'H': 2, 'O': 1}],
[True, True, False, False],
[2.0, 32895.0, 25188.0, 25422.0]],
"PG5 combustion"
),
(
[[{'C': 597000000, 'H': 744000002, 'N': 225000000, 'O': 390000001, 'P': 60000000}, {'O': 2}, {'C': 1, 'O': 2}, {'H': 2, 'O': 1}, {'N': 1, 'O': 2}, {'P': 4, 'O': 10}],
[True, True, False, False, False, False],
[1.0, 888000000.0, 597000000.0, 372000001.0, 225000000.0, 15000000.0]],
"Y chromosome DNA combustion: C597000000H744000002N225000000O390000001P60000000 + O2 = CO2 + H2O + NO2 + P4O10"
),
(
[[{'C': 597000000, 'H': 744000002, 'N': 225000000, 'O': 390000001, 'P': 60000000}, {'H': 2, 'O': 2}, {'C': 1, 'O': 2}, {'H': 2, 'O': 1}, {'N': 2}, {'P': 2, 'O': 5}],
[True, True, False, False, False, False],
[1.0, 1326000000.0, 597000000.0, 1698000000.0, 112500000.0, 30000000.0]],
"Y chromosome with H2O2: C597000000H744000002N225000000O390000001P60000000 + H2O2 = CO2 + H2O + N2 + P2O5"
),
(
[[{'C': 597000000, 'H': 744000002, 'N': 225000000, 'O': 390000001, 'P': 60000000}, {'H': 2, 'S': 1, 'O': 4}, {'C': 1, 'O': 2}, {'H': 2, 'O': 1}, {'N': 1, 'O': 2}, {'H': 1, 'P': 1, 'O': 4}, {'S': 1, 'O': 2}],
[True, True, False, False, False, False, False],
[1.0, 1836000000.0, 597000000.0, 2178000000.0, 225000000.0, 60000000.0, 1836000000.0]],
"Y chromosome with H2SO4: C597000000H744000002N225000000O390000001P60000000 + H2SO4 = CO2 + H2O + NO2 + HPO4 + SO2"
),
(
[[{'C': 597000000, 'H': 744000002, 'N': 225000000, 'O': 390000001, 'P': 60000000}, {'O': 2}, {'C': 1, 'O': 1}, {'H': 2, 'O': 1}, {'N': 1, 'O': 2}, {'P': 4, 'O': 10}],
[True, True, False, False, False, False],
[1.0, 589500000.0, 597000000.0, 372000001.0, 225000000.0, 15000000.0]],
"Y chromosome incomplete combustion: C597000000H744000002N225000000O390000001P60000000 + O2 = CO + H2O + NO2 + P4O10"
),
(
[[{'C': 597000000, 'H': 744000002, 'N': 225000000, 'O': 390000001, 'P': 60000000}, {'Cl': 2}, {'C': 1, 'Cl': 4}, {'H': 1, 'Cl': 1}, {'N': 1, 'Cl': 3}, {'P': 1, 'Cl': 5}, {'O': 2}],
[True, True, False, False, False, False, False],
[2.0, 4107000002.0, 1194000000.0, 1488000004.0, 450000000.0, 120000000.0, 390000001.0]],
"Y chromosome chlorination: C597000000H744000002N225000000O390000001P60000000 + Cl2 = CCl4 + HCl + NCl3 + PCl5 + O2"
)
]
@pytest.mark.parametrize("test_case,test_name", ill_conditioned_stoich_test_cases)
def test_balance_stoichiometry_ill_conditioned(test_case, test_name):
"""Test stoichiometry balancing for ill-conditioned reactions."""
atomss, statuses, products = test_case
matrix = stoichiometric_matrix(atomss, statuses)
products_calc = balance_stoichiometry(matrix, rounding=11)
assert_close1d(products_calc, products)
atol = sum(abs(v) for r in matrix for v in r)*1e-10
check_reaction_balance(matrix, products_calc, atol=atol)

stoich_test_cases = [
[[{'Hg': 1, 'O': 1}, {'Hg': 1}, {'O': 2}], [True, False, False], [2.0, 2.0, 1.0]],
[[{'Cl': 2}, {'C': 3, 'H': 6}, {'C': 3, 'Cl': 1, 'H': 5}, {'Cl': 1, 'H': 1}],
[True, True, False, False, False],
Expand Down Expand Up @@ -1719,13 +1794,19 @@ def test_balance_stoichiometry():
[True, True, False, False],
[1.0, 1.0, 1.0, 2.0]],

]
]

for atomss, statuses, products in test_cases:
@pytest.mark.parametrize("test_case", stoich_test_cases)
def test_balance_stoichiometry(test_case):
atomss, statuses, products = test_case
for settings in [{'rounding': 9, 'allow_fractional': False},
{'rounding': 16, 'allow_fractional': True}]:
matrix = stoichiometric_matrix(atomss, statuses)
products_calc = balance_stoichiometry(matrix)
check_reaction_balance(matrix, products_calc)
assert_close1d(products_calc, products)
products_calc = balance_stoichiometry(matrix, **settings)
if not settings['allow_fractional']:
# when we allow fractions we stil have valid ratios but they do not match the hardcoded answers
assert_close1d(products_calc, products)
check_reaction_balance(matrix, products_calc, atol=1e-12)


def test_round_to_significant():
Expand Down

0 comments on commit 6e6d3b4

Please sign in to comment.