Commit d0def804 authored by eladnoor's avatar eladnoor

small changes regarding pint

parent d51b5b54
Pipeline #46348133 canceled with stage
in 36 seconds
......@@ -25,7 +25,7 @@
# THE SOFTWARE.
from equilibrator_api import ComponentContribution, Reaction, ureg, kJ_per_mol
from equilibrator_api import ComponentContribution, Reaction, Q_
import logging
import argparse
import pandas as pd
......@@ -45,19 +45,35 @@ if __name__ == '__main__':
parser.add_argument(
'outfile', type=argparse.FileType('w'),
help='path to output file')
parser.add_argument('--i', type=float,
help='ionic strength in M',
default=0.1)
parser.add_argument('--i', type=str,
help='ionic strength',
default="0.1 M")
parser.add_argument('--t', type=str,
help='ionic strength',
default="298.15 K")
parser.add_argument('--ph', type=float, help='pH level', default=7.0)
logging.getLogger().setLevel(logging.WARNING)
args = parser.parse_args()
p_h = args.ph
ionic_strength = args.i * ureg.molar
p_h = Q_(args.ph)
assert p_h.check(None)
sys.stderr.write('pH = %.1f\n' % p_h)
ionic_strength = Q_(args.i)
if ionic_strength.check(None):
ionic_strength *= Q_("M")
else:
assert ionic_strength.check("[concentration]")
temperature = Q_(args.t)
if temperature.check(None):
temperature *= Q_("K")
else:
assert temperature.check("[temperature]")
sys.stderr.write('pH = %s\n' % p_h)
sys.stderr.write('I = %s\n' % ionic_strength)
sys.stderr.write('T = %s\n' % temperature)
infile_lines = list(filter(None, map(str.strip, args.infile.readlines())))
......@@ -69,20 +85,20 @@ if __name__ == '__main__':
else:
reactions = list(map(Reaction.parse_formula, infile_lines))
eq_api = ComponentContribution(p_h=p_h, ionic_strength=ionic_strength)
eq_api = ComponentContribution(p_h=p_h, ionic_strength=ionic_strength,
temperature=temperature)
result_df = pd.DataFrame(data=list(zip(infile_lines, reactions)),
columns=["reaction (original)",
"reaction (parsed)"])
result_df['pH'] = p_h
result_df['I [M]'] = float(ionic_strength/ureg.molar)
dG0_prime, U = eq_api.standard_dg_prime_multi(reactions)
result_df["ΔG'° [kJ/mol]"] = [float(dg0/kJ_per_mol) for dg0 in dG0_prime]
result_df["σ(ΔG'°) [kJ/mol]"] = [
float(np.sqrt(u)/kJ_per_mol) for u in U.diagonal().flat]
result_df["ln(Reversibility index)"] = [
float(eq_api.ln_reversibility_index(r)) for r in reactions]
result_df['I'] = ionic_strength
standard_dg_prime, U = eq_api.standard_dg_prime_multi(reactions)
result_df["ΔG'°"] = list(standard_dg_prime.flat)
result_df["σ(ΔG'°)"] = list(map(np.sqrt, U.diagonal().flat))
result_df["ln(Reversibility index)"] = list(map(
eq_api.ln_reversibility_index, reactions))
result_df["comment"] = ""
non_balanced = result_df["reaction (parsed)"].apply(
......@@ -91,6 +107,5 @@ if __name__ == '__main__':
result_df.loc[non_balanced, "ln(Reversibility index)"] = np.nan
result_df.loc[non_balanced, "comment"] = "reaction is not chemically balanced"
result_df.to_csv(args.outfile)
args.outfile.flush()
......@@ -29,7 +29,7 @@ from ._version import get_versions
__version__ = get_versions()['version']
del get_versions
from equilibrator_cache import Compound, ureg, kJ_per_mol, R
from equilibrator_cache import Compound, ureg, Q_, R
from equilibrator_cache.compound_cache import ccache
from equilibrator_api.bounds import Bounds
from equilibrator_api.pathway import Pathway
......
......@@ -110,7 +110,7 @@ class BaseBounds(object):
b2log = lambda b: log(b / standard_concentration)
return map(b2log, lbs), map(b2log, ubs)
@ureg.check(None, None, 'molar', 'molar')
@ureg.check(None, None, '[concentration]', '[concentration]')
def SetBounds(self, compound: Compound, lb: float, ub: float):
"""Set bounds for a specific key
......
......@@ -41,7 +41,7 @@ class ComponentContribution(object):
holds default conditions for compounds in the different phases,
"""
@ureg.check(None, None, None, 'molar', 'kelvin')
@ureg.check(None, None, None, '[concentration]', '[temperature]')
def __init__(
self,
p_h: float = default_pH,
......
......@@ -27,11 +27,11 @@
from equilibrator_cache.thermodynamic_constants import (
FARADAY, POSSIBLE_REACTION_ARROWS, R, Rlog10, default_I, default_pH,
default_T, kJ_per_mol, physiological_concentration, standard_concentration,
ureg)
default_T, physiological_concentration, standard_concentration,
Q_)
default_pMg = 14.0
default_pMg = Q_("14.0")
default_phase = 'aqueous'
default_conc_lb = 1e-6 * ureg.molar
default_conc_ub = 1e-2 * ureg.molar
default_conc_lb = Q_("1e-6 M")
default_conc_ub = Q_("1e-2 M")
......@@ -34,7 +34,7 @@ import pandas as pd
from optlang.glpk_interface import Constraint, Model, Objective, Variable
from .settings import (
R, default_T, kJ_per_mol, physiological_concentration,
R, default_T, Q_, physiological_concentration,
standard_concentration)
......@@ -175,10 +175,10 @@ class PathwayMDFData(object):
@property
def reaction_plot(self) -> plt.Figure:
cumulative_dgms = np.cumsum([0] + self.reaction_df.dGm.tolist())
cumulative_dgms *= float(R * default_T / kJ_per_mol)
cumulative_dgms *= float(R * default_T / Q_("kJ/mol"))
cumulative_dgs = np.cumsum([0] + self.reaction_df.dGr.tolist())
cumulative_dgs *= float(R * default_T / kJ_per_mol)
cumulative_dgs *= float(R * default_T / Q_("kJ/mol"))
xticks = np.arange(0, self.reaction_df.shape[0] + 1)
xticklabels = [''] + self.reaction_df.reaction_id.tolist()
......
......@@ -30,13 +30,17 @@ import warnings
import pytest
from equilibrator_api import (
ComponentContribution, Reaction, ccache, kJ_per_mol, ureg)
ComponentContribution, Reaction, ccache, Q_)
def assert_with_unit(obs_dg: float, exp_dg: float, abs: float = 1e-1,
unit: str = "kJ/mol"):
assert obs_dg.check(unit)
assert float(obs_dg / Q_(unit)) == pytest.approx(float(exp_dg), abs=abs)
@pytest.fixture(scope="module")
def comp_contribution() -> ComponentContribution:
return ComponentContribution(p_h = 7.0,
ionic_strength = 0.1 * ureg.molar)
return ComponentContribution(p_h = Q_("7"), ionic_strength = Q_("0.1M"))
def test_gibbs_energy_atp_hydrolysis(comp_contribution):
......@@ -45,22 +49,22 @@ def test_gibbs_energy_atp_hydrolysis(comp_contribution):
# ATP + H2O = ADP + Pi
formula = "KEGG:C00002 + KEGG:C00001 = KEGG:C00008 + KEGG:C00009"
reaction = Reaction.parse_formula(formula)
compound_to_conc = {ccache.get_compound('KEGG:C00002'): 1e-3 * ureg.molar,
ccache.get_compound('KEGG:C00008'): 1e-2 * ureg.molar,
ccache.get_compound('KEGG:C00009'): 1e-2 * ureg.molar}
compound_to_conc = {ccache.get_compound('KEGG:C00002'): Q_("1mM"),
ccache.get_compound('KEGG:C00008'): Q_("10mM"),
ccache.get_compound('KEGG:C00009'): Q_("10mM")}
standard_dg_prime, dg_uncertainty = comp_contribution.standard_dg_prime(
reaction)
assert (standard_dg_prime / (-26.4 * kJ_per_mol)) == pytest.approx(1.0, abs=1e-2)
assert (dg_uncertainty / (0.3 * kJ_per_mol)) == pytest.approx(1.0,
abs=1e-1)
assert_with_unit(standard_dg_prime, -26.4)
assert_with_unit(dg_uncertainty, 0.3)
dg_prime, _ = comp_contribution.dg_prime(reaction, compound_to_conc)
assert (dg_prime / (-32.1 * kJ_per_mol)) == pytest.approx(1.0, abs=1e-2)
assert_with_unit(dg_prime, -32.1)
physiological_dg_prime, _ = comp_contribution.physiological_dg_prime(reaction)
assert (physiological_dg_prime / (-43.5 * kJ_per_mol)) == pytest.approx(1.0, abs=1e-2)
physiological_dg_prime, _ = comp_contribution.physiological_dg_prime(
reaction)
assert_with_unit(physiological_dg_prime, -43.5)
def test_gibbs_energy_pyruvate_decarboxylase(comp_contribution):
......@@ -70,10 +74,11 @@ def test_gibbs_energy_pyruvate_decarboxylase(comp_contribution):
formula = "KEGG:C00022 = KEGG:C00084 + KEGG:C00011"
reaction = Reaction.parse_formula(formula)
dG0_prime, dG0_uncertainty = comp_contribution.standard_dg_prime(reaction)
standard_dg_prime, dg_uncertainty = comp_contribution.standard_dg_prime(
reaction)
assert dG0_prime / (-18.1 * kJ_per_mol) == pytest.approx(1.0, abs=1e-1)
assert dG0_uncertainty / (3.3 * kJ_per_mol) == pytest.approx(1.0, abs=1e-1)
assert_with_unit(standard_dg_prime, -18.4)
assert_with_unit(dg_uncertainty, 3.3)
def test_reduction_potential(comp_contribution):
......@@ -83,11 +88,18 @@ def test_reduction_potential(comp_contribution):
formula = "KEGG:C00036 = KEGG:C00149"
reaction = Reaction.parse_formula(formula)
E0_prime, E0_uncertainty = comp_contribution.standard_e_prime(reaction)
standard_e_prime, e_uncertainty = comp_contribution.standard_e_prime(
reaction)
assert_with_unit(standard_e_prime, -175.0, unit="mV")
assert_with_unit(e_uncertainty, 3.3, unit="mV")
def test_unresolved_reactions(comp_contribution):
formula = "KEGG:C09844 + KEGG:C00003 + KEGG:C00001 => KEGG:C03092 + KEGG:C00004"
reaction = Reaction.parse_formula(formula)
_, dg_uncertainty = comp_contribution.standard_dg_prime(reaction)
expected_E0_prime = -178e-3 * ureg.volts
expected_E0_uncertainty = 3.3e-3 * ureg.volts
assert float(dg_uncertainty / Q_("kJ/mol")) > 1e4
assert E0_prime / expected_E0_prime == pytest.approx(1.0, abs=1e-2)
assert E0_uncertainty / expected_E0_uncertainty == pytest.approx(1.0,
abs=1e-2)
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment