Commit 841b7f41 authored by eladnoor's avatar eladnoor

adding units using pint

parent feb9e463
......@@ -112,17 +112,21 @@ Then inside the python3 shell:
from equilibrator_api import Pathway
pp = Pathway.from_sbtab('example.tsv')
mdf_res = pp.calc_mdf()
print('MDF = ', mdf_res.mdf)
mdf_res.mdf_plot.show()
print('MDF = %.2f RT', mdf_res.mdf)
mdf_res.reaction_plot.show()
```
## Dependencies:
- python >= 3.6
- equilibrator-cache
- component-contribution
- sbtab
- numpy
- scipy
- optlang
- pandas
- nltk
- pyparsing
- sbtab
- matplotlib
- quilt
- pint
......@@ -63,9 +63,9 @@ if __name__ == '__main__':
reactions.append(Reaction({}))
continue
equilibrator = ComponentContribution(pH=args.ph, ionic_strength=args.i)
cc = ComponentContribution(pH=args.ph, ionic_strength=args.i)
dG0_prime, U = equilibrator.dG0_prime_multi(reactions)
dG0_prime, U = cc.standard_dg_prime_multi(reactions)
writer = csv.writer(args.outfile)
header = ['reaction', 'pH', 'ionic strength [M]', 'dG\'0 [kJ/mol]',
......@@ -77,7 +77,7 @@ if __name__ == '__main__':
if r.is_empty():
row += [nan, nan, nan, 'reaction is empty']
elif r.check_full_reaction_balancing():
ln_RI = r.calculate_reversibility_index_from_dG0_prime(dg0)
ln_RI = cc.reversibility_index(r)
row += ['%.2f' % dg0, '%.2f' % sqrt(u), '%.2f' % ln_RI, '']
else:
row += [nan, nan, nan, 'reaction is not chemically balanced']
......
......@@ -50,8 +50,8 @@ if __name__ == '__main__':
output_pdf = PdfPages(args.outfile)
mdf_res = pp.calc_mdf()
output_pdf.savefig(mdf_res.conc_plot)
output_pdf.savefig(mdf_res.mdf_plot)
output_pdf.savefig(mdf_res.compound_plot)
output_pdf.savefig(mdf_res.reaction_plot)
output_pdf.close()
rxn_df = mdf_res.reaction_df
......
......@@ -25,7 +25,6 @@ keywords =
[options]
zip_safe = True
install_requires =
importlib_resources
numpy>=1.15.2
scipy>=1.1.0
optlang>=1.4.3
......@@ -34,6 +33,7 @@ install_requires =
pyparsing>=2.2.0
sbtab>=0.9.49
matplotlib>=3.0.0
pint>=0.9
requests
equilibrator-cache>=0.1.4
component-contribution>=0.1.0
......
......@@ -30,12 +30,12 @@ from typing import Dict, Iterable, TextIO, Tuple
import pandas as pd
import pkg_resources
from equilibrator_cache import Compound
from equilibrator_cache import Compound, ureg
from equilibrator_cache.compound_cache import ccache
from numpy import log
from equilibrator_api import settings
from equilibrator_api.concs import ConcentrationConverter
from equilibrator_api.settings import (
default_conc_lb, default_conc_ub, standard_concentration)
class InvalidBounds(Exception):
......@@ -49,10 +49,6 @@ class BaseBounds(object):
"""Returns a (deep) copy of self"""
raise NotImplementedError
def GetRange(self):
"""Returns a 2-tuple of the concentration range"""
return None
def GetLowerBound(self, compound: Compound):
"""Get the lower bound for this key
......@@ -112,8 +108,10 @@ class BaseBounds(object):
items are Numpy arrays of dimensions 1xlen(keys)
"""
lbs, ubs = self.GetBounds(compounds)
return map(log, lbs), map(log, ubs)
b2log = lambda b: log(b / standard_concentration)
return map(b2log, lbs), map(b2log, ubs)
@ureg.check(None, None, 'molar', 'molar')
def SetBounds(self, compound: Compound, lb: float, ub: float):
"""Set bounds for a specific key
......@@ -133,12 +131,12 @@ class Bounds(BaseBounds):
Contains upper and lower bounds for various keys
Allows for defaults
"""
@ureg.check(None, None, None, '[concentration]', '[concentration]')
def __init__(self,
lower_bounds: Dict[Compound, float] = None,
upper_bounds: Dict[Compound, float] = None,
default_lb: float = settings.DEFAULT_CONC_LB,
default_ub: float = settings.DEFAULT_CONC_UB) -> object:
default_lb: float = default_conc_lb,
default_ub: float = default_conc_ub) -> object:
"""Initialize the bounds object.
Args:
......@@ -149,17 +147,22 @@ class Bounds(BaseBounds):
"""
self.lower_bounds = lower_bounds or dict()
self.upper_bounds = upper_bounds or dict()
for b in self.lower_bounds.values():
assert b.check('[concentration]')
for b in self.upper_bounds.values():
assert b.check('[concentration]')
self.default_lb = default_lb
self.default_ub = default_ub
self.c_range = (self.default_lb, self.default_ub)
@classmethod
@ureg.check(None, None, '[concentration]', '[concentration]')
def from_csv(cls,
f: TextIO,
default_lb: float = settings.DEFAULT_CONC_LB,
default_ub: float = settings.DEFAULT_CONC_UB) -> object:
"""Read bounds from a .csv file
default_lb: float = default_conc_lb,
default_ub: float = default_conc_ub) -> object:
"""Read bounds from a .csv file. Assume that all the concentrations
are given in units of molar.
Args:
f: an open .csv file stream
......@@ -172,19 +175,21 @@ class Bounds(BaseBounds):
for row in bounds_df.itertuples():
compound = ccache.get_compound(row.compound_id)
lbs[compound] = row.lb
ubs[compound] = row.ub
lbs[compound] = row.lb * ureg.molar
ubs[compound] = row.ub * ureg.molar
bounds = Bounds(lbs, ubs, default_lb, default_ub)
bounds._check_bounds()
return bounds
@classmethod
@ureg.check(None, None, '[concentration]', '[concentration]',
'[concentration]')
def from_dataframe(cls,
df: pd.DataFrame,
unit: str = 'Molar',
default_lb: float = settings.DEFAULT_CONC_LB,
default_ub: float = settings.DEFAULT_CONC_UB) -> object:
concentration_unit: ureg.Unit,
default_lb: float = default_conc_lb,
default_ub: float = default_conc_ub) -> object:
"""Read bounds from a pandas.DataFrame
Args:
......@@ -199,13 +204,8 @@ class Bounds(BaseBounds):
row = df.loc[idx]
cid = row['Compound:Identifiers:kegg.compound']
compound = ccache.get_compound_from_registry('kegg', cid)
ub = float(row['Concentration:Max'])
lb = float(row['Concentration:Min'])
ub = ConcentrationConverter.to_molar_string(ub, unit)
lb = ConcentrationConverter.to_molar_string(lb, unit)
ubs[compound] = ub
lbs[compound] = lb
lbs[compound] = ureg(row['Concentration:Min']) * concentration_unit
ubs[compound] = ureg(row['Concentration:Max']) * concentration_unit
bounds = Bounds(lbs, ubs, default_lb, default_ub)
bounds._check_bounds()
......@@ -220,19 +220,19 @@ class Bounds(BaseBounds):
def _check_bounds(self) -> None:
if self.default_lb > self.default_ub:
msg = (
'Default lower bound %.2g > default upper bound %.2g' %
(self.default_lb, self.default_ub))
raise InvalidBounds(msg)
raise InvalidBounds(
'Default lower bound %s > default upper bound %s' %
(self.default_lb, self.default_ub)
)
for compound in self.upper_bounds:
lb = self.GetLowerBound(compound)
ub = self.GetUpperBound(compound)
if lb > ub:
msg = (
'Invalid bounds for %s: lower bound %f > upper bound %f' %
(compound, lb, ub))
raise InvalidBounds(msg)
raise InvalidBounds(
'Invalid bounds for %s: lower bound %s > upper bound %s' %
(compound, lb, ub)
)
def Copy(self) -> object:
"""
......@@ -244,18 +244,13 @@ class Bounds(BaseBounds):
self.default_lb,
self.default_ub)
def GetRange(self):
"""Returns a 2-tuple of the concentration range"""
return self.c_range
def GetLowerBound(self, key: Compound):
"""Get the lower bound for this key
Args:
key: a compound
"""
val = self.lower_bounds.get(key) or self.default_lb
return val
return self.lower_bounds.get(key, self.default_lb)
def GetUpperBound(self, key: Compound):
"""Get the upper bound for this key.
......@@ -263,15 +258,18 @@ class Bounds(BaseBounds):
Args:
key: a compound
"""
val = self.upper_bounds.get(key) or self.default_ub
return val
COFACTORS_FNAME = pkg_resources.resource_filename(
'equilibrator_api', 'data/cofactors.csv')
with open(COFACTORS_FNAME, 'r') as fp:
DEFAULT_BOUNDS = Bounds.from_csv(fp)
#df = DEFAULT_BOUNDS.to_dataframe()
#df.to_csv('cofactors.csv')
return self.upper_bounds.get(key, self.default_ub)
# the default is generated only at the first call of GetDefaultBounds()
# the reason to do this, is because looking up all the cofactor compounds
# in the cache takes about 20 seconds
DEFAULT_BOUNDS = None
@staticmethod
def GetDefaultBounds():
if Bounds.DEFAULT_BOUNDS is None:
COFACTORS_FNAME = pkg_resources.resource_filename(
'equilibrator_api', 'data/cofactors.csv')
with open(COFACTORS_FNAME, 'r') as fp:
Bounds.DEFAULT_BOUNDS = Bounds.from_csv(fp)
return Bounds.DEFAULT_BOUNDS
......@@ -29,7 +29,7 @@ from typing import Dict, List, Tuple
import numpy as np
from component_contribution import GibbsEnergyPredictor
from equilibrator_cache import Compound
from equilibrator_cache import Compound, ureg
from equilibrator_api import settings
from equilibrator_api.phased_reaction import PhasedReaction
......@@ -41,26 +41,30 @@ class ComponentContribution(object):
holds default conditions for compounds in the different phases,
"""
@ureg.check(None, None, None, 'molar', 'kelvin')
def __init__(
self,
pH: float = settings.DEFAULT_PH,
pMg: float = settings.DEFAULT_PMG,
ionic_strength: float = settings.DEFAULT_IONIC_STRENGTH
self,
p_h: float = settings.default_pH,
p_mg: float = settings.default_pMg,
ionic_strength: float = settings.default_I,
temperature: float = settings.default_T
) -> object:
self.predictor = GibbsEnergyPredictor()
self.pH = pH
self.p_h = p_h
self.ionic_strength = ionic_strength
self.pMg = pMg
self.T = settings.DEFAULT_TEMP
self.p_mg = p_mg
self.temperature = temperature
@property
def RT(self) -> float:
return settings.R * self.T
return settings.R * self.temperature
def dG0_prime_multi(self,
reactions: List[PhasedReaction]
) -> Tuple[np.array, np.array]:
def standard_dg_prime_multi(
self,
reactions: List[PhasedReaction],
raise_exception: bool = False
) -> Tuple[np.array, np.array]:
"""Calculate the transformed reaction energies of a list of reactions.
:return: a tuple of two arrays. the first is a 1D NumPy array
......@@ -71,15 +75,19 @@ class ComponentContribution(object):
U as the covariance of a Gaussian used for sampling
(where dG0_cc is the mean of that Gaussian).
"""
return self.predictor.get_dG0_r_prime_multi(
return self.predictor.standard_dg_prime_multi(
reactions,
pH=self.pH,
p_h=self.p_h,
ionic_strength=self.ionic_strength,
T=self.T,
raise_exception=False
T=self.temperature,
raise_exception=raise_exception
)
def dG0_prime(self, reaction: PhasedReaction) -> Tuple[float, float]:
def standard_dg_prime(
self,
reaction: PhasedReaction,
raise_exception: bool = False
) -> Tuple[float, float]:
"""Calculate the transformed reaction energies of a reaction.
:param reaction: the input Reaction object
......@@ -87,18 +95,19 @@ class ComponentContribution(object):
calculate the confidence interval, use the range -1.96 to 1.96 times
this value
"""
return self.predictor.get_dG0_r_prime(
return self.predictor.standard_dg_prime(
reaction,
pH=self.pH,
p_h=self.p_h,
ionic_strength=self.ionic_strength,
T=self.T,
raise_exception=False
temperature=self.temperature,
raise_exception=raise_exception
)
def dG_prime(
self,
reaction: PhasedReaction,
compound_to_conc: Dict[Compound, float]
def dg_prime(
self,
reaction: PhasedReaction,
compound_to_conc: Dict[Compound, float],
raise_exception: bool = False
) -> Tuple[float, float]:
"""
Calculate the dG'0 of a single reaction
......@@ -114,12 +123,17 @@ class ComponentContribution(object):
(which is the value shown on eQuilibrator)
"""
dG0_r_prime, dG0_uncertainty = self.dG0_prime(reaction)
dG_r_prime = dG0_r_prime + (self.RT *
reaction.dG_correction(compound_to_conc))
return dG_r_prime, dG0_uncertainty
dg_prime, dg_uncertainty = self.standard_dg_prime(
reaction,
raise_exception=raise_exception
)
dg_prime += self.RT * reaction.dg_correction(compound_to_conc)
return dg_prime, dg_uncertainty
def dGm_prime(self, reaction: PhasedReaction) -> Tuple[float, float]:
def physiological_dg_prime(
self,
reaction: PhasedReaction
) -> Tuple[float, float]:
"""
Calculate the dG'm of a single reaction (i.e. if all reactants are
at 1 mM, except H2O)
......@@ -128,15 +142,17 @@ class ComponentContribution(object):
reaction - an object of type Reaction
Returns:
dGm_r_prime - estimated Gibbs free energy of reaction
dGm_uncertainty - standard deviation of estimation, multiply
by 1.96 to get a 95% confidence interval
(which is the value shown on eQuilibrator)
physiological_dg_prime - estimated Gibbs free energy of reaction
dg_uncertainty - standard deviation of estimation, multiply
by 1.96 to get a 95% confidence interval
(which is the value shown on eQuilibrator)
"""
dG0_r_prime, dG0_uncertainty = self.dG0_prime(reaction)
dGm_r_prime = dG0_r_prime + self.RT * reaction.dGm_correction()
return dGm_r_prime, dG0_uncertainty
physiological_dg_prime, dg_uncertainty = self.standard_dg_prime(
reaction)
physiological_dg_prime += self.RT * \
reaction.physiological_dg_correction()
return physiological_dg_prime, dg_uncertainty
def reversibility_index(self, reaction: PhasedReaction) -> float:
"""
......@@ -148,15 +164,15 @@ class ComponentContribution(object):
by 1.96 to get a 95% confidence interval
(which is the value shown on eQuilibrator)
"""
dGm_r_prime, dG0_uncertainty = self.dGm_prime(reaction)
physiological_dg_prime, _ = self.physiological_dg_prime(reaction)
abs_sum_coeff = reaction._GetAbsSumCoeff()
if abs_sum_coeff == 0:
return np.inf
ln_RI = (2.0 / abs_sum_coeff) * dGm_r_prime / self.RT
ln_RI = (2.0 / abs_sum_coeff) * physiological_dg_prime / self.RT
return ln_RI
def E0_prime(self, reaction: PhasedReaction) -> float:
def standard_e_prime(self, reaction: PhasedReaction) -> Tuple[float, float]:
"""
Calculate the E'0 of a single half-reaction
"""
......@@ -167,16 +183,16 @@ class ComponentContribution(object):
raise ValueError('this is not a half-reaction, '
'electrons are balanced')
dG0_prime, dG0_uncertainty = self.dG0_prime(reaction)
dG0_prime, dG0_uncertainty = self.standard_dg_prime(reaction)
E0_prime_mV = 1000.0 * -dG0_prime / (n_e*settings.FARADAY)
E0_uncertainty = 1000.0 * dG0_uncertainty / (n_e*settings.FARADAY)
E0_prime = -dG0_prime / (n_e*settings.FARADAY) # in Volts
E0_uncertainty = dG0_uncertainty / (n_e*settings.FARADAY) # in Volts
return E0_prime_mV, E0_uncertainty
return E0_prime, E0_uncertainty
def dG0_analysis(self, reaction: PhasedReaction) -> List[Dict[str, object]]:
return self.predictor.get_dG0_r_analysis(reaction,
raise_exception=False)
return self.predictor.get_dg_analysis(reaction,
raise_exception=False)
def IsUsingGroupContributions(self, reaction: PhasedReaction) -> bool:
return self.predictor.is_using_group_contribution(reaction,
......
......@@ -30,15 +30,15 @@ import logging
from typing import Callable, Iterable, List, TextIO, Tuple
import numpy as np
from equilibrator_cache import ureg
from equilibrator_cache.compound_cache import Compound, ccache
from equilibrator_cache.reaction import (
create_stoichiometric_matrix_from_reactions)
from sbtab import SBtab
from equilibrator_api.bounds import DEFAULT_BOUNDS, Bounds
from equilibrator_api.bounds import Bounds
from equilibrator_api.component_contribution import ComponentContribution
from equilibrator_api.phased_reaction import PhasedReaction
from equilibrator_api.settings import RT
from equilibrator_api.thermo_models import PathwayThermoModel
......@@ -67,8 +67,8 @@ class Pathway(object):
def __init__(self,
reactions: List[PhasedReaction],
fluxes: Iterable[float],
dG0_r_primes: Iterable[float],
dG0_r_uncertainty: np.array = None,
standard_dg_primes: Iterable[float],
dg_uncertainties: np.array = None,
bounds: Bounds = None
) -> object:
"""Initialize.
......@@ -76,8 +76,8 @@ class Pathway(object):
Args:
reactions: a list of gibbs.reaction.Reaction objects.
fluxes: numpy.array of relative fluxes in same order as reactions.
dG0_r_primes: reaction energies.
dG0_r_uncertainty: covariance matrix of uncertainties in dG'0
standard_dg_primes: reaction energies (in units of RT)
dg_uncertainties: covariance matrix of uncertainties (in RT units)
bounds: bounds on metabolite concentrations.
Uses default bounds if None provided.
"""
......@@ -85,20 +85,22 @@ class Pathway(object):
Nr = len(reactions)
self.fluxes = np.array(list(fluxes))
self.dG0_r_prime = np.array(list(dG0_r_primes))
self.standard_dg_primes = np.array(list(standard_dg_primes))
assert self.fluxes.shape == (Nr, )
assert self.dG0_r_prime.shape == (Nr, )
assert self.standard_dg_primes.shape == (Nr,)
if dG0_r_uncertainty is None:
self.dG0_r_uncertainty = np.zeros((len(reactions), len(reactions)))
if dg_uncertainties is None:
self.dg_uncertainties = np.zeros((len(reactions), len(reactions)))
else:
assert dG0_r_uncertainty.shape == (Nr, Nr)
self.dG0_r_uncertainty = dG0_r_uncertainty
assert dg_uncertainties.shape == (Nr, Nr)
self.dg_uncertainties = dg_uncertainties
phys_corr = np.array(
[float(r.physiological_dg_correction()) for r in reactions])
self.physiological_dg_prime = self.standard_dg_primes + phys_corr
self.dGm_r_prime = (self.dG0_r_prime +
np.array([r.dGm_correction() for r in reactions]))
if bounds is None:
self._bounds = DEFAULT_BOUNDS.Copy()
self._bounds = Bounds.GetDefaultBounds().Copy()
else:
self._bounds = bounds.Copy()
......@@ -115,7 +117,7 @@ class Pathway(object):
S_T = self.S.T.values
S_inv = np.linalg.pinv(S_T)
null_proj = np.eye(self.S.shape[1]) - S_T @ S_inv
projected = null_proj * self.dG0_r_prime.T
projected = null_proj * self.standard_dg_primes.T
if (projected > 1e-8).any():
raise ViolatesFirstLaw(
'Supplied reaction dG values are inconsistent '
......@@ -169,12 +171,11 @@ class Pathway(object):
reactions.append(rxn)
fluxes.append(flux)
equilibrator = ComponentContribution(pH=pH,
ionic_strength=ionic_strength)
dG0_r_primes, dG0_uncertainties = zip(*map(equilibrator.dG0_prime,
cc = ComponentContribution(pH=pH, ionic_strength=ionic_strength)
standard_dg_primes, _ = zip(*map(cc.standard_dg_prime,
reactions))
dG0_r_primes = list(dG0_r_primes)
return Pathway(reactions, fluxes, dG0_r_primes, bounds=bounds)
standard_dg_primes = list(standard_dg_primes)
return Pathway(reactions, fluxes, standard_dg_primes, bounds=bounds)
@staticmethod
def get_compounds(reactions: List[PhasedReaction]) -> List[Compound]:
......@@ -235,9 +236,12 @@ class Pathway(object):
compound_to_name = bounds_df.Compound.to_dict()
name_to_compound = dict(map(reversed, compound_to_name.items()))
bounds_unit = sbtabdoc.get_sbtab_by_id(
bounds_unit_str = sbtabdoc.get_sbtab_by_id(
'ConcentrationConstraint').get_attribute('Unit')
bounds_unit = ureg(bounds_unit_str)
assert bounds_unit.check('[concentration]'), "Bounds must be in " \
"concentration units, " \
f"not '{bounds_unit_str}'"
bounds = Bounds.from_dataframe(bounds_df, bounds_unit)
reactions = []
......@@ -259,14 +263,14 @@ class Pathway(object):
# grab rows containing keqs.
keqs = keqs_df[keqs_df['QuantityType'] == 'equilibrium constant']
reaction_keqs = dict(zip(keqs['Reaction'], keqs['Value']))
dgs = [-RT * np.log(float(reaction_keqs[rid]))
for rid in reaction_ids]
# Manually set the delta G values on the reaction objects
for dg, rxn in zip(dgs, reactions):
rxn._dg0_prime = dg
# calculate the standard Gibbs reaction energies from the equilibrium
# constants, in units of RT
standard_dg_primes = [-np.log(float(reaction_keqs[rid]))
for rid in reaction_ids]
pp = Pathway(reactions, fluxes_ordered, dgs, bounds=bounds)
pp = Pathway(reactions, fluxes_ordered, standard_dg_primes,
bounds=bounds)
# override the compound names with the ones in the SBtab
pp.set_compound_names(compound_to_name.get)
......
......@@ -31,11 +31,12 @@ import re
import numpy
import pyparsing
from .settings import (
POSSIBLE_COMPOUND_BODYCHARS, POSSIBLE_COMPOUND_INITCHARS,
POSSIBLE_REACTION_ARROWS)
from .settings import POSSIBLE_REACTION_ARROWS
POSSIBLE_COMPOUND_INITCHARS = pyparsing.alphanums + "()"
POSSIBLE_COMPOUND_BODYCHARS = pyparsing.alphanums + "-+,()[]'_"
class ParseError(Exception):
pass
......
......@@ -25,30 +25,13 @@
# THE SOFTWARE.
from numpy import log
from pyparsing import alphanums
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)
POSSIBLE_REACTION_ARROWS = ('<->', '<=>', '-->', '<--', # 3-character arrows
'=>', '<=', '->', '<-', # 2-character arrows
'=', '⇌', '⇀', '⇋', '↽') # 1-character arrows
POSSIBLE_COMPOUND_INITCHARS = alphanums + "()"
POSSIBLE_COMPOUND_BODYCHARS = alphanums + "-+,()[]'_"
R = 8.31e-3 # kJ/(K*mol)
DEFAULT_TEMP = 298.15 # K
DEFAULT_IONIC_STRENGTH = 0.1 # mM
DEFAULT_PH = 7.0
DEFAULT_PMG = 14.0
DEFAULT_PHASE = 'aqueous'
RT = R * DEFAULT_TEMP
RTlog10 = RT * log(10)
DEBYE_HUECKLE_A = 2.91482
DEBYE_HUECKLE_B = 1.6
MG_FORMATION_ENERGY = -455.3 # kJ/mol, formation energy of Mg2+
FARADAY = 96.485 # kC/mol
DEFAULT_CONC_LB = 1e-6
DEFAULT_CONC_UB = 1e-2
default_pMg = 14.0
default_phase = 'aqueous'
default_conc_lb = 1e-6 * ureg.molar
default_conc_ub = 1e-2 * ureg.molar
This diff is collapsed.
......@@ -30,6 +30,7 @@ import os
import warnings
import pytest
from equilibrator_cache import kJ_per_mol, ureg
from equilibrator_api import PhasedReaction
from equilibrator_api.pathway import Pathway
......@@ -53,7 +54,7 @@ def test_mdf():
"2 KEGG:C00002 + 2 KEGG:C00001 + 2 KEGG:C00469 + 2 KEGG:C00011")
assert net_rxn == reference_net_rxn
assert pp_mdf_data.mdf == pytest.approx(1.69, abs=1e-2)
assert pp_mdf_data.mdf == pytest.approx(0.682, rel=1e-2)
shadow_prices = pp_mdf_data.reaction_df.set_index(
'reaction_id').shadow_price
......@@ -66,12 +67,12 @@ def test_mdf():
compound_prices = pp_mdf_data.compound_df.set_index('compound').shadow_price
assert compound_prices['h2o'] == pytest.approx(0.0, abs=1e-2)
assert compound_prices['nad'] == pytest.approx(1.24, abs=1e-2)
assert compound_prices['nadh'] == pytest.approx(-1.24, abs=1e-2)
assert compound_prices['nad'] == pytest.approx(0.5, abs=1e-2)
assert compound_prices['nadh'] == pytest.approx(-0.5, abs=1e-2)
# test the mdf_plot and conc_plot run without crashing
# test that the plotting functions run without crashing
try:
pp_mdf_data.mdf_plot
pp_mdf_data.conc_plot
pp_mdf_data.reaction_plot
pp_mdf_data.compound_plot
except:
pytest.fail("MDF plot function raises exception")
......@@ -30,61 +30,68 @@ import os
import warnings
import pytest
from equilibrator_cache import kJ_per_mol, ureg
from equilibrator_cache.compound_cache import ccache
from equilibrator_api.component_contribution import ComponentContribution
from equilibrator_api.phased_reaction import PhasedReaction
TEST_DIR = os.path.dirname(
os.path.abspath(inspect.getfile(inspect.currentframe())))
@pytest.fixture(scope="module")
def comp_contribution() -> ComponentContribution:
return ComponentContribution(p_h = 7.0,
ionic_strength = 0.1 * ureg.molar)