Commit 247fd884 authored by eladnoor's avatar eladnoor

flake8

parent 911597ee
......@@ -40,11 +40,12 @@ if __name__ == '__main__':
parser = argparse.ArgumentParser(
description='Calculate potentials for a number of reactions.')
parser.add_argument('infile', type=argparse.FileType('r'),
help='path to input file containing reactions')
help='path to input file containing reactions')
parser.add_argument('outfile', type=argparse.FileType('w'),
help='path to output file')
help='path to output file')
parser.add_argument('--plaintext', action='store_true',
help='indicate that reactions are given in plain text (not accessions)')
help="'indicate that reactions are given in plain text"
" (not accessions)")
parser.add_argument('--ph', type=str,
help='pH level',
default="7.0")
......@@ -80,11 +81,8 @@ if __name__ == '__main__':
infile_lines = list(filter(None, map(str.strip, args.infile.readlines())))
if args.plaintext:
raise Exception("ReactionMatcher needs to be reimplemented")
#reaction_matcher = ReactionMatcher()
#reactions = list(map(reaction_matcher.match, infile_lines))
raise NotImplementedError("ReactionMatcher needs to be reimplemented")
else:
reactions = list(map(Reaction.parse_formula, infile_lines))
......@@ -105,10 +103,11 @@ if __name__ == '__main__':
result_df["comment"] = ""
non_balanced = result_df["reaction (parsed)"].apply(
lambda r: r.is_balanced(raise_exception=False) == False
lambda r: not r.is_balanced(raise_exception=False)
)
result_df.loc[non_balanced, "ln(Reversibility index)"] = np.nan
result_df.loc[non_balanced, "comment"] = "reaction is not chemically balanced"
result_df.loc[non_balanced, "comment"] = ("reaction is not chemically "
"balanced")
result_df.to_csv(args.outfile)
args.outfile.flush()
......@@ -28,7 +28,6 @@
import argparse
import logging
import pandas as pd
from matplotlib.backends.backend_pdf import PdfPages
from equilibrator_api import Pathway
......@@ -56,7 +55,6 @@ if __name__ == '__main__':
pdf.savefig(mdf_res.compound_plot)
pdf.savefig(mdf_res.reaction_plot)
net_rxn = pp.net_reaction()
rxn_df = mdf_res.reaction_df
cpd_df = mdf_res.compound_df
......@@ -62,6 +62,7 @@ universal = 1
max-line-length = 80
exclude =
__init__.py
_version.py
[pydocstyle]
match_dir = equilibrator_api
......
......@@ -30,9 +30,12 @@ del get_versions
from equilibrator_cache import Compound, ureg, Q_, R
from equilibrator_cache.compound_cache import ccache
from equilibrator_cache.thermodynamic_constants import (
FARADAY, POSSIBLE_REACTION_ARROWS, Q_, R, default_I, default_pH,
default_pMg, default_T, physiological_concentration,
standard_concentration, ureg)
from equilibrator_api.bounds import Bounds
from equilibrator_api.pathway import Pathway
from equilibrator_api.phased_reaction import PhasedReaction as Reaction
from equilibrator_api.component_contribution import ComponentContribution
from equilibrator_api.thermo_models import PathwayThermoModel
......@@ -31,10 +31,9 @@ import pandas as pd
import pkg_resources
from numpy import log
from equilibrator_api.settings import (
default_conc_lb, default_conc_ub, standard_concentration)
from .settings import default_conc_lb, default_conc_ub
from . import Q_, Compound, ccache, ureg
from . import Q_, Compound, ccache, ureg, standard_concentration
class InvalidBounds(Exception):
......@@ -121,6 +120,10 @@ class BaseBounds(object):
lbs, ubs = zip(*bounds)
return lbs, ubs
@staticmethod
def b2log(b):
log(b / standard_concentration)
def GetLnBounds(
self,
compounds: Iterable[Compound]
......@@ -135,8 +138,8 @@ class BaseBounds(object):
items are Numpy arrays of dimensions 1xlen(keys)
"""
lbs, ubs = self.GetBounds(compounds)
b2log = lambda b: log(b / standard_concentration)
return map(b2log, lbs), map(b2log, ubs)
return map(self.b2log, lbs), map(self.b2log, ubs)
def GetLnLowerBounds(
self,
......@@ -151,8 +154,7 @@ class BaseBounds(object):
An iterable of the log-scaled lower bounds
"""
lbs = self.GetLowerBounds(compounds)
b2log = lambda b: log(b / standard_concentration)
return map(b2log, lbs)
return map(self.b2log, lbs)
def GetLnUpperBounds(
self,
......@@ -167,9 +169,7 @@ class BaseBounds(object):
An iterable of the log-scaled upper bounds
"""
ubs = self.GetUpperBounds(compounds)
b2log = lambda b: log(b / standard_concentration)
return map(b2log, ubs)
return map(self.b2log, ubs)
@ureg.check(None, None, '[concentration]', '[concentration]')
def SetBounds(self, compound: Compound, lb: float, ub: float):
......@@ -223,7 +223,7 @@ class Bounds(BaseBounds):
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
default_lb: default lower bound to if no specific one is provided
......
......@@ -29,11 +29,10 @@ from typing import Dict, List, Tuple
import numpy as np
from component_contribution import GibbsEnergyPredictor
from equilibrator_cache.reaction import Reaction
from . import Compound, ureg
from .phased_reaction import AQUEOUS_PHASE, PHASED_COMPOUND_DICT, PhasedReaction
from .settings import FARADAY, R, default_I, default_pH, default_pMg, default_T
from .phased_reaction import PhasedReaction
from . import FARADAY, R, default_I, default_pH, default_pMg, default_T
class ComponentContribution(object):
......@@ -75,12 +74,10 @@ class ComponentContribution(object):
U as the covariance of a Gaussian used for sampling
(where dG0_cc is the mean of that Gaussian).
"""
rxn_aquoues, dg_prime_nonaq = zip(*
[reaction.split_aqueous_phase(
p_h=self.p_h,
ionic_strength=self.ionic_strength,
temperature=self.temperature)
for reaction in reactions])
rxn_dg_pairs = map(lambda r: r.split_aqueous_phase(
p_h=self.p_h, ionic_strength=self.ionic_strength,
temperature=self.temperature), reactions)
rxn_aquoues, dg_prime_nonaq = zip(*rxn_dg_pairs)
dg_prime_nonaq = np.array(dg_prime_nonaq)
standard_dg_prime, dg_sigma = self.predictor.standard_dg_prime_multi(
......@@ -160,8 +157,8 @@ class ComponentContribution(object):
physiological_dg_prime, dg_uncertainty = self.standard_dg_prime(
reaction)
physiological_dg_prime += self.RT * \
reaction.physiological_dg_correction()
physiological_dg_prime += \
self.RT * reaction.physiological_dg_correction()
return physiological_dg_prime, dg_uncertainty
def ln_reversibility_index(self, reaction: PhasedReaction) -> float:
......@@ -192,8 +189,8 @@ class ComponentContribution(object):
dG0_prime, dG0_uncertainty = self.standard_dg_prime(reaction)
E0_prime = -dG0_prime / (n_e*FARADAY) # in Volts
E0_uncertainty = dG0_uncertainty / (n_e*FARADAY) # in Volts
E0_prime = -dG0_prime / (n_e*FARADAY) # in Volts
E0_uncertainty = dG0_uncertainty / (n_e*FARADAY) # in Volts
return E0_prime, E0_uncertainty
......
......@@ -41,8 +41,9 @@ from . import Compound
from .bounds import Bounds
from .component_contribution import ComponentContribution
from .phased_reaction import PhasedReaction
from .settings import (
Q_, R, default_I, default_pH, default_pMg, default_T, strip_units)
from . import (
Q_, R, default_I, default_pH, default_pMg, default_T)
from .settings import strip_units
from .thermo_models import PathwayThermoModel
......@@ -179,10 +180,10 @@ class Pathway(object):
mapping: Callable[[Compound], str]
) -> None:
"""Use alternative compound names for outputs such as plots
Args:
mapping: a dictionary mapping compounds to their names in the model
"""
self.compound_names = list(map(mapping, self.S.index))
......@@ -264,24 +265,28 @@ class Pathway(object):
return pp
@classmethod
def from_sbtab(self, sbtab) -> object:
"""
read the sbtab file (can be a filename or file handel)
and use it to initialize the Pathway
"""
def open_sbtab(cls, sbtab) -> SBtab.SBtabDocument:
if type(sbtab) == SBtab.SBtabDocument:
sbtabdoc = sbtab
return sbtab
elif type(sbtab) == str:
sbtabdoc = SBtab.read_csv(sbtab, "pathway")
return SBtab.read_csv(sbtab, "pathway")
elif isinstance(sbtab, IOBase):
sbtab_contents = sbtab.read()
if type(sbtab_contents) == bytes:
sbtab_contents = sbtab_contents.decode('utf-8')
sbtabdoc = SBtab.SBtabDocument("pathway", sbtab_contents,
"unnamed_sbtab.tsv")
else:
raise ValueError("sbtab must be either a file name or a Stream "
"object")
return SBtab.SBtabDocument("pathway", sbtab_contents,
"unnamed_sbtab.tsv")
raise ValueError("sbtab must be either a file name or a Stream "
"object")
@classmethod
def from_sbtab(cls, sbtab) -> object:
"""
read the sbtab file (can be a filename or file handel)
and use it to initialize the Pathway
"""
sbtabdoc = cls.open_sbtab(sbtab)
table_ids = ['ConcentrationConstraint', 'Reaction', 'RelativeFlux',
'Parameter']
dfs = []
......
......@@ -42,12 +42,13 @@ SOLID_PHASE = 's'
LIQUID_PHASE = 'l'
PhaseInfo = namedtuple("phase_info", "name "
"standard_abundance "
"physiolical_abundance "
"dimensionality")
"standard_abundance "
"physiolical_abundance "
"dimensionality")
PHASE_INFO_DICT = {
AQUEOUS_PHASE: PhaseInfo("aqueous",Q_("1 M"),Q_("1 mM"),"[concentration]"),
AQUEOUS_PHASE: PhaseInfo("aqueous", Q_("1 M"), Q_("1 mM"),
"[concentration]"),
GAS_PHASE: PhaseInfo("gas", Q_("1 bar"), Q_("1 mbar"), "[pressure]"),
SOLID_PHASE: PhaseInfo("solid", None, None, "[dimensionless]"),
LIQUID_PHASE: PhaseInfo("liquid", None, None, "[dimensionless]")
......@@ -73,17 +74,17 @@ MicroSpecie = namedtuple("microspecie", "standard_dgf num_protons charge")
# None, use the value from component_contribution
PHASED_COMPOUND_DICT = {
("MNXM4", GAS_PHASE): MicroSpecie(
standard_dgf=Q_("0 kJ/mol"), num_protons=0, charge=0), # O2
standard_dgf=Q_("0 kJ/mol"), num_protons=0, charge=0), # O2
("MNXM195", GAS_PHASE): MicroSpecie(
standard_dgf=Q_("0 kJ/mol"), num_protons=2, charge=0), # H2
standard_dgf=Q_("0 kJ/mol"), num_protons=2, charge=0), # H2
("MNXM724", GAS_PHASE): MicroSpecie(
standard_dgf=Q_("0 kJ/mol"), num_protons=0, charge=0), # N2
standard_dgf=Q_("0 kJ/mol"), num_protons=0, charge=0), # N2
("MNXM13", GAS_PHASE): MicroSpecie(
standard_dgf=Q_("-394.36 kJ/mol"), num_protons=0, charge=0), # CO2
standard_dgf=Q_("-394.36 kJ/mol"), num_protons=0, charge=0), # CO2
("MNXM10881", GAS_PHASE): MicroSpecie(
standard_dgf=Q_("-137.17 kJ/mol"), num_protons=0, charge=0), # CO
standard_dgf=Q_("-137.17 kJ/mol"), num_protons=0, charge=0), # CO
("MNXM89", SOLID_PHASE): MicroSpecie(
standard_dgf=Q_("0 kJ/mol"), num_protons=0, charge=0), # S8
standard_dgf=Q_("0 kJ/mol"), num_protons=0, charge=0), # S8
}
......@@ -219,11 +220,11 @@ class PhasedReaction(Reaction):
Takes an unbalanced reaction and converts it into an oxidation
reaction, by adding H2O, O2, Pi, CO2, and NH4+ to both sides.
"""
balancing_ids = ['MNX:MNXM2', # H2O
'MNX:MNXM4', # O2
'MNX:MNXM9', # Pi
'MNX:MNXM13', # CO2
'MNX:MNXM15', # NH4+
balancing_ids = ['MNX:MNXM2', # H2O
'MNX:MNXM4', # O2
'MNX:MNXM9', # Pi
'MNX:MNXM13', # CO2
'MNX:MNXM15', # NH4+
]
S = ccache.get_element_data_frame(balancing_ids)
......@@ -261,6 +262,7 @@ class PhasedReaction(Reaction):
"""
return PhasedReaction({kegg_id: -1}).balance_by_oxidation()
if __name__ == '__main__':
import sys
......
# The MIT License (MIT)
#
# Copyright (c) 2013 Weizmann Institute of Science
# Copyright (c) 2018 Institute for Molecular Systems Biology,
# ETH Zurich
# Copyright (c) 2018 Novo Nordisk Foundation Center for Biosustainability,
# Technical University of Denmark
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
import logging
import re
import numpy
import pyparsing
from .settings import POSSIBLE_REACTION_ARROWS
POSSIBLE_COMPOUND_INITCHARS = pyparsing.alphanums + "()"
POSSIBLE_COMPOUND_BODYCHARS = pyparsing.alphanums + "-+,()[]'_"
class ParseError(Exception):
pass
def _parsedCompound(c_list):
"""Always put a stoichiometric coefficient with a compound."""
if len(c_list) == 2:
return c_list[0], c_list[1]
return 1, c_list[0]
def _MakeReactionSideParser():
"""Builds a parser for a side of a reaction."""
# Coefficients are usually integral, but they can be floats or fractions too.
int_coeff = pyparsing.Word(pyparsing.nums)
float_coeff = pyparsing.Word(pyparsing.nums + '.' + pyparsing.nums)
frac_coeff = int_coeff + '/' + int_coeff
int_coeff.setParseAction(lambda i:int(i[0]))
float_coeff.setParseAction(lambda t:float(t[0]))
frac_coeff.setParseAction(lambda f:float(f[0])/float(f[2]))
coeff = pyparsing.Or([int_coeff, float_coeff, frac_coeff])
optional_coeff = pyparsing.Optional(coeff)
compound_separator = pyparsing.Literal('+').suppress()
compound_name_component = pyparsing.Word(
initChars=POSSIBLE_COMPOUND_INITCHARS,
bodyChars=POSSIBLE_COMPOUND_BODYCHARS)
compound_name = pyparsing.Forward()
compound_name << (compound_name_component + pyparsing.ZeroOrMore(compound_name_component))
compound_name.setParseAction(lambda s: ' '.join(s))
compound_with_coeff = pyparsing.Forward()
compound_with_coeff << ((optional_coeff + compound_name) | compound_name)
compound_with_coeff.setParseAction(_parsedCompound)
compound_with_coeff.setResultsName("compound")
compound_with_separator = pyparsing.Forward()
compound_with_separator << (compound_with_coeff + compound_separator)
reaction_side = pyparsing.Forward()
reaction_side << (pyparsing.ZeroOrMore(compound_with_separator) +
compound_with_coeff)
reaction_side.setParseAction(lambda l: [l])
reaction_side.setResultsName("reaction_side")
return reaction_side
def _MakeReactionParser():
"""Builds a pyparsing-based recursive descent parser for chemical reactions."""
reaction_side = _MakeReactionSideParser()
side_separators = [pyparsing.Literal(s) for s in POSSIBLE_REACTION_ARROWS]
side_separator = pyparsing.Or(side_separators).suppress()
reaction = pyparsing.Forward()
reaction << (reaction_side + side_separator + reaction_side)
return reaction
class ParsedReactionQuery(object):
"""A parsed reaction query."""
def __init__(self, substrates=None, products=None):
"""Initialize the ParsedReaction object.
Args:
reactants: a list of tuples for the reactants.
products: a list of tuples for the products.
"""
self.substrates = substrates or []
self.products = products or []
def __eq__(self, other):
"""Equality test."""
r = frozenset(self.substrates)
p = frozenset(self.products)
o_r = frozenset(other.substrates)
o_p = frozenset(other.products)
reactants_diff = r.symmetric_difference(o_r)
products_diff = p.symmetric_difference(o_p)
if not reactants_diff and not products_diff:
return True
return False
def __str__(self):
joined_rs = ['%s %s' % (numpy.abs(c),r) for c,r in self.substrates]
joined_ps = ['%s %s' % (numpy.abs(c),p) for c,p in self.products]
return '%s => %s' % (' + '.join(joined_rs), ' + '.join(joined_ps))
class QueryParser(object):
"""Parses search queries."""
REACTION_PATTERN = u'.*(' + '|'.join(POSSIBLE_REACTION_ARROWS) + ').*'
REACTION_MATCHER = re.compile(REACTION_PATTERN)
def __init__(self):
"""Initialize the parser."""
self._rparser = _MakeReactionParser()
def is_reaction_query(self, query):
"""Returns True if this query is likely to be a reaction query.
Args:
query: the query string.
"""
m = self.REACTION_MATCHER.match(query.strip())
return m is not None
def parse_reaction_query(self, query):
"""Parse the query as a reaction.
Args:
query: the query string.
Returns:
An initialized ParsedReaction object, or None if parsing failed.
"""
try:
results = self._rparser.parseString(query)
substrates, products = results
logging.debug('substrates = %s' % str(substrates))
logging.debug('products = %s' % str(products))
return ParsedReactionQuery(substrates, products)
except pyparsing.ParseException as msg:
logging.error('Failed to parse query %s', query)
raise ParseError(msg)
if __name__ == '__main__':
qp = QueryParser()
# example of a standard reaction (ATP hydrolysis)
parsed_query = qp.parse_reaction_query('ATP + H2O ↽ ADP + orthophosphate')
print(parsed_query.substrates)
This diff is collapsed.
......@@ -25,15 +25,14 @@
# THE SOFTWARE.
import numpy as np
from equilibrator_cache.thermodynamic_constants import (
FARADAY, POSSIBLE_REACTION_ARROWS, Q_, R, default_I, default_pH, default_T,
physiological_concentration, standard_concentration, ureg)
from . import Q_
default_pMg = Q_("14.0")
default_phase = 'aqueous'
default_conc_lb = Q_("1e-6 M")
default_conc_ub = Q_("1e-2 M")
def strip_units(v: np.array) -> np.array:
return np.array(list(map(float, v.flat))).reshape(v.shape)
......@@ -34,7 +34,7 @@ import pandas as pd
from matplotlib.lines import Line2D
from optlang.glpk_interface import Constraint, Model, Objective, Variable
from .settings import (
from . import (
Q_, R, default_T, physiological_concentration, standard_concentration, ureg)
......@@ -45,7 +45,7 @@ class PathwayMDFData(object):
pathway,
B: float,
I_dir: np.array,
l: np.array,
lnC: np.array,
y: np.array,
reaction_prices: np.array,
compound_prices: np.array
......@@ -54,7 +54,7 @@ class PathwayMDFData(object):
self.mdf = B * R * default_T
self.bounds = pathway._bounds.Copy()
concentrations = np.exp(l) * standard_concentration
concentrations = np.exp(lnC) * standard_concentration
standard_dg_primes = pathway.standard_dg_primes
physiological_dg_prime = PathwayMDFData.calc_physiological_dg(pathway)
......@@ -64,9 +64,10 @@ class PathwayMDFData(object):
optimized_dg_prime += pathway.dg_sigma @ y
# adjust dG to flux directions and convert to kJ/mol
standard_dg_primes = (I_dir @ standard_dg_primes) * R * default_T
physiological_dg_prime = (I_dir @ physiological_dg_prime) * R * default_T
optimized_dg_prime = (I_dir @ optimized_dg_prime) * R * default_T
RT = R * pathway.comp_contrib.temperature
standard_dg_primes = (I_dir @ standard_dg_primes) * RT
physiological_dg_prime = (I_dir @ physiological_dg_prime) * RT
optimized_dg_prime = (I_dir @ optimized_dg_prime) * RT
# all dG values are in units of RT, so we convert them to kJ/mol
reaction_data = zip(pathway.reaction_ids(),
......@@ -152,7 +153,8 @@ class PathwayMDFData(object):
data_df.loc[data_df.shadow_price != 0, 'color'] = 'red'
# a sub-DataFrame of only the metabolites with fixed concentrations
data_df.loc[data_df.lower_bound == data_df.upper_bound, 'color'] = 'grey'
data_df.loc[
data_df.lower_bound == data_df.upper_bound, 'color'] = 'grey'
compound_fig, ax = plt.subplots(1, 1, figsize=(9, 6))
ax.xaxis.set_units(ureg.molar)
......@@ -200,7 +202,7 @@ class PathwayMDFData(object):
cumulative_dgms = np.array([float(x/Q_('kJ/mol')) for x in
cumulative_dgms]) * Q_('kJ/mol')
cumulative_dgs = np.array([float(x/Q_('kJ/mol')) for x in
cumulative_dgs]) * Q_('kJ/mol')
cumulative_dgs]) * Q_('kJ/mol')
xticks = 0.5 + np.arange(self.reaction_df.shape[0])
xticklabels = self.reaction_df.reaction_id.tolist()
......@@ -244,9 +246,9 @@ class PathwayThermoModel(object):
# Make sure dG0_r' is the right size
assert pathway.standard_dg_primes.shape == (self.Nr, ), \
'standard dG required for all reactions'
'standard dG required for all reactions'
assert pathway.dg_sigma.shape == (self.Nr, self.Nr), \
'uncertainty in dG required for all reactions'
'uncertainty in dG required for all reactions'
self.I_dir = np.diag(np.sign(self.pathway.fluxes))
self.Nr_active = sum(self.pathway.fluxes.T != 0)
......@@ -314,11 +316,11 @@ class PathwayThermoModel(object):
b1 = -self.I_dir[inds] @ self.pathway.standard_dg_primes
b2 = np.ones(self.Nr) * self.stdev_factor
A = np.vstack([np.hstack([ A11, A12, A13]), # driving force
np.hstack([ A21, A22, A23]), # covariance var ub
np.hstack([-A21, A22, A23]), # covariance var lb
np.hstack([ A31, A32, A33]), # log conc ub
np.hstack([ A31, -A32, A33])]) # log conc lb
A = np.vstack([np.hstack([A11, A12, A13]), # driving force
np.hstack([A21, A22, A23]), # covariance var ub
np.hstack([-A21, A22, A23]), # covariance var lb
np.hstack([A31, A32, A33]), # log conc ub
np.hstack([A31, -A32, A33])]) # log conc lb
b = np.hstack([b1,
b2,
......@@ -353,9 +355,9 @@ class PathwayThermoModel(object):
y = [Variable("y%d" % i) for i in range(self.Nr)]
# ln-concentration variables
l = [Variable("l%d" % i) for i in range(self.Nc)]
lnC = [Variable("l%d" % i) for i in range(self.Nc)]
return A, b, c, y, l
return A, b, c, y, lnC
def _GetDualVariablesAndConstants(
self
......@@ -380,19 +382,21 @@ class PathwayThermoModel(object):
Does not set the objective function... leaves that to the caller.
Returns:
the linear problem object, and the three types of variables as arrays
the linear problem object, and the three types of variables as
arrays
"""
A, b, c, y, l = self._GetPrimalVariablesAndConstants()
A, b, c, y, lnC = self._GetPrimalVariablesAndConstants()
B = Variable('mdf')
x = y + l + [B]
x = y + lnC + [B]
lp = Model(name="MDF_PRIMAL")
cnstr_names = ["driving_force_%02d" % j for j in range(self.Nr_active)] + \
["covariance_var_ub_%02d" % j for j in range(self.Nr)] + \
["covariance_var_lb_%02d" % j for j in range(self.Nr)] + \
["log_conc_ub_%02d" % j for j in range(self.Nc)] + \
["log_conc_lb_%02d" % j for j in range(self.Nc)]
cnstr_names = \
[f"driving_force_{j:02d}" for j in range(self.Nr_active)] + \
[f"covariance_var_ub_{j:02d}" for j in range(self.Nr)] + \
[f"covariance_var_lb_{j:02d}" for j in range(self.Nr)] + \
[f"log_conc_ub_{j:02d}" for j in range(self.Nc)] + \
[f"log_conc_lb_{j:02d}" for j in range(self.Nc)]
constraints = []
for j in range(A.shape[0]):
......@@ -405,7 +409,7 @@ class PathwayThermoModel(object):
row = [c[i] * x[i] for i in range(c.shape[0])]
lp.objective = Objective(sum(row), direction='max')
return lp, y, l, B
return lp, y, lnC, B
def _MakeMDFProblemDual(
self
......@@ -426,13 +430,13 @@ class PathwayThermoModel(object):
cnstr_names = ["y_%02d" % j for j in range(self.Nr)] + \
["l_%02d" % j for j in range(self.Nc)] + \
["MDF"]
constraints = []
for i in range(A.shape[1]):
row = [A[j, i] * x[j] for j in range(A.shape[0])]
constraints.append(Constraint(sum(row), lb=c[i], ub=c[i],
name=cnstr_names[i]))
lp.add(constraints)
row = [b[i] * x[i] for i in range(A.shape[0])]
......@@ -453,15 +457,15 @@ class PathwayThermoModel(object):
"""
def get_primal_array(l):
return np.array([v.primal for v in l], ndmin=1)
lp_primal, y, l, B = self._MakeMDFProblem()
lp_primal, y, lnC, B = self._MakeMDFProblem()
if lp_primal.optimize() != 'optimal':
logging.warning('LP status %s', lp_primal.status)
raise Exception("Cannot solve MDF primal optimization problem")
y = get_primal_array(y) # covariance eigenvalue prefactors
l = get_primal_array(l) # log concentrations
y = get_primal_array(y) # covariance eigenvalue prefactors
lnC = get_primal_array(lnC) # log concentrations
B = lp_primal.variables['mdf'].primal
lp_dual, w, g, z, u = self._MakeMDFProblemDual()
......@@ -471,8 +475,8 @@ class PathwayThermoModel(object):
primal_obj = lp_primal.objective.value
dual_obj = lp_dual.objective.value
if abs(primal_obj - dual_obj) > 1e-3:
raise Exception("Primal != Dual (%.5f != %.5f)"
% (primal_obj, dual_obj))
raise Exception(f"Primal != Dual ({primal_obj:.5f} != "
f"{dual_obj:.5f}")