Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fixing some bugs and active development #28

Open
wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Fixing some bugs and active development #28

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
13 changes: 7 additions & 6 deletions setup.cfg
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,15 +36,16 @@ keywords =
[options]
zip_safe = True
install_requires =
cobra
cobra==0.23
depinfo
optlang<1.4.6
numpy
optlang>1.5
numpy==1.23
scipy
pandas
equilibrator-api==0.3.2b7
component-contribution==0.3.2b4
equilibrator-cache==0.3.2b2
equilibrator-api
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It would be nice to pin it, since the different equilibrator packages usually need to coordinated. This version set works for me with this branch:

equilibrator-api==0.6.0
equilibrator-cache==0.6.0
component-contribution==0.6.0

I think stating only equilibrator-api>=0.6.0 and letting pip do the dependency resolution might also be fine.

component-contribution
equilibrator-cache
python-libsbml

python_requires = >=3.6
tests_require =
Expand Down
18 changes: 14 additions & 4 deletions src/multitfa/analysis/variability.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,13 @@

import numpy as np
from pandas import DataFrame, Series
import os

cwd = os.getcwd()

logs_dir = cwd + os.sep + "logs"
if not os.path.exists(logs_dir):
os.makedirs(logs_dir)


def variability(model_variability, variable_list=None):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have been using this branch with CPLEX for a couple of days, including the variability, and it's been working great so far.

Expand Down Expand Up @@ -217,6 +224,9 @@ def variability_legacy_cplex(
if not os.path.exists(tmp_dir):
os.makedirs(tmp_dir)

# Log file for debugging
cplexlog = open(logs_dir + os.sep + "cplex.log", "a")

# Instantiate Cplex model
cplex_model = Cplex()
rand_str = "".join(choices(string.ascii_lowercase + string.digits, k=6))
Expand All @@ -227,10 +237,10 @@ def variability_legacy_cplex(
model.cplex_interface.write(temp_filename)
cplex_model.read(temp_filename)

cplex_model.set_log_stream(None)
cplex_model.set_error_stream(None)
cplex_model.set_warning_stream(None)
cplex_model.set_results_stream(None)
cplex_model.set_log_stream(cplexlog)
cplex_model.set_error_stream(cplexlog)
cplex_model.set_warning_stream(cplexlog)
cplex_model.set_results_stream(cplexlog)

if params:
# print("lol")
Expand Down
16 changes: 8 additions & 8 deletions src/multitfa/core/compound.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,7 @@ def concentration_variable(self):
optlang.interface.variable
optlang interface variable for concentration of metabolite. if metabolite is not linked to model, then None
"""
if self.model is not None:
if self.model:
conc_var = "lnc_{}".format(self.id)
return self.model.variables[conc_var]
else:
Expand All @@ -63,7 +63,7 @@ def concentration_min(self):
try:
return self._concentration_min
except AttributeError:
if self.model.compartment_info is not None:
if not self.model.compartment_info.empty:
try:
self._concentration_min = float(
self.model.compartment_info["c_min"][self.compartment]
Expand All @@ -82,13 +82,13 @@ def concentration_min(self, value):
lb=np.log(value), ub=np.log(self.concentration_max)
)
# Update the Gurobi and Cplex interface bounds
if self.model.gurobi_interface is not None:
if self.model.gurobi_interface:
self.model.gurobi_interface.getVarByName(
self.concentration_variable.name
).LB = np.log(value)
self.model.gurobi_interface.update()

if self.model.cplex_interface is not None:
if self.model.cplex_interface:
self.model.cplex_interface.variables.set_lower_bounds(
self.concentration_variable.name, np.log(value)
)
Expand All @@ -105,7 +105,7 @@ def concentration_max(self):
try:
return self._concentration_max
except AttributeError:
if self.model.compartment_info is not None:
if not self.model.compartment_info.empty:
try:
self._concentration_max = float(
self.model.compartment_info["c_max"][self.compartment]
Expand All @@ -124,13 +124,13 @@ def concentration_max(self, value):
lb=np.log(self.concentration_min), ub=np.log(value)
)
# Update the Gurobi and cplex interface bounds
if self.model.gurobi_interface is not None:
if self.model.gurobi_interface:
self.model.gurobi_interface.getVarByName(
self.concentration_variable.name
).UB = np.log(value)
self.model.gurobi_interface.update()

if self.model.cplex_interface is not None:
if self.model.cplex_interface:
self.model.cplex_interface.variables.set_upper_bounds(
self.concentration_variable.name, np.log(value)
)
Expand Down Expand Up @@ -166,7 +166,7 @@ def delG_err_variable(self):
optlang.interface.variable
metabolite Gibbs energy error variable, if metabolite is not associated with model, then None
"""
if self.model is not None:
if self.model:
err_var = "dG_err_{}".format(self.id)
return self.model.variables[err_var]
else:
Expand Down
77 changes: 66 additions & 11 deletions src/multitfa/core/reaction.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@
from cobra import Reaction
from numpy import zeros
from six import iteritems

from math import isinf
from ..util.thermo_constants import FARADAY


Expand Down Expand Up @@ -53,7 +53,7 @@ def delG_forward(self):
An optlang variable for the Gibbs energy of the forward half reaction or None if the reaction is not associated with the model or reaction is in the list of thermodynamic excluded reactions.
"""
var_name = "dG_{}".format(self.forward_variable.name)
if self.model is not None:
if self.model:
if self.id not in self.model.Exclude_reactions:
return self.model.variables[var_name]
else:
Expand All @@ -71,7 +71,7 @@ def delG_reverse(self):
An optlang variable for the Gibbs energy of the reverse half reaction or None if the reaction is not associated with the model or reaction is in the list of thermodynamic excluded reactions.
"""
var_name = "dG_{}".format(self.reverse_variable.name)
if self.model is not None:
if self.model:
if self.id not in self.model.Exclude_reactions:
return self.model.variables[var_name]
else:
Expand All @@ -89,7 +89,7 @@ def indicator_forward(self):
An optlang binary variable of forward half reaction or None if the reaction is not associated with the model or eaction is in the list of thermodynamic excluded reactions.
"""
var_name = "indicator_{}".format(self.forward_variable.name)
if self.model is not None:
if self.model:
if self.id not in self.model.Exclude_reactions:
return self.model.variables[var_name]
else:
Expand All @@ -107,7 +107,7 @@ def indicator_reverse(self):
An optlang binary variable of reverse half reaction or None if the reaction is not associated with the model or eaction is in the list of thermodynamic excluded reactions.
"""
var_name = "indicator_{}".format(self.reverse_variable.name)
if self.model is not None:
if self.model:
if self.id not in self.model.Exclude_reactions:
return self.model.variables[var_name]
else:
Expand All @@ -129,7 +129,7 @@ def directionality_constraint(self):
forward_cons = "directionality_{}".format(self.forward_variable.name)
reverse_cons = "directionality_{}".format(self.reverse_variable.name)

if self.model is not None:
if self.model:
if self.id not in self.model.Exclude_reactions:
return (
self.model.constraints[forward_cons],
Expand All @@ -154,7 +154,7 @@ def indicator_constraint(self):
forward_cons = "ind_{}".format(self.forward_variable.name)
reverse_cons = "ind_{}".format(self.reverse_variable.name)

if self.model is not None:
if self.model:
if self.id not in self.model.Exclude_reactions:
return (
self.model.constraints[forward_cons],
Expand All @@ -181,7 +181,7 @@ def delG_constraint(self):
forward_cons = "delG_{}".format(self.forward_variable.name)
reverse_cons = "delG_{}".format(self.reverse_variable.name)

if self.model is not None:
if self.model:
if self.id not in self.model.Exclude_reactions:
return (
self.model.constraints[forward_cons],
Expand Down Expand Up @@ -223,7 +223,8 @@ def delG_transport(self):
try:
return self._delG_transport
except AttributeError:
self._delG_transport = self.calculate_delG_transport()
electro_delG, proton_delG = self.calculate_delG_transport()
self._delG_transport = electro_delG + proton_delG
return self._delG_transport

@delG_transport.setter
Expand Down Expand Up @@ -261,6 +262,60 @@ def calculate_transport_charge(self):
{key: sum(val) for key, val in n_proton.items()},
)

def update_variable_bounds(self):
"""Updating the method from cobra to update the flux bounds in cplex_interface. Currently not doing it for gurobi_interface as I don't have a licence to test this"""
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have a Gurobi license, could give it a try and make a PR to this branch or in a follow-up after this is merged. Or would you rather stop support for gurobi? My use case for gurobi so far has been that the constraint relaxation is more flexible, but I have managed to get CPLEX to do what I want. Not sure if the "Debugging" part of the docs would still work, although I don't think that gurobi's computeIIS function takes the quadratic constraint into account anyways.

if self.model is None:
return
super().update_variable_bounds()
if self._lower_bound > 0:
if self.model.cplex_interface:
self.model.cplex_interface.variables.set_lower_bounds(
self.forward_variable.name,
None if isinf(self._lower_bound) else self._lower_bound,
)
self.model.cplex_interface.variables.set_upper_bounds(
self.forward_variable.name,
None if isinf(self._upper_bound) else self._upper_bound,
)

self.model.cplex_interface.variables.set_lower_bounds(
self.reverse_variable.name, 0
)
self.model.cplex_interface.variables.set_upper_bounds(
self.reverse_variable.name, 0
)
elif self._upper_bound < 0:
self.model.cplex_interface.variables.set_lower_bounds(
self.forward_variable.name, 0
)
self.model.cplex_interface.variables.set_upper_bounds(
self.forward_variable.name, 0
)

self.model.cplex_interface.variables.set_lower_bounds(
self.reverse_variable.name,
None if isinf(self._upper_bound) else -self._upper_bound,
)
self.model.cplex_interface.variables.set_upper_bounds(
self.reverse_variable.name,
None if isinf(self._lower_bound) else -self._lower_bound,
)
else:
self.model.cplex_interface.variables.set_lower_bounds(
self.forward_variable.name, 0
)
self.model.cplex_interface.variables.set_upper_bounds(
self.forward_variable.name,
None if isinf(self._upper_bound) else self._upper_bound,
)
self.model.cplex_interface.variables.set_lower_bounds(
self.reverse_variable.name, 0
)
self.model.cplex_interface.variables.set_upper_bounds(
self.reverse_variable.name,
None if isinf(self._lower_bound) else -self._lower_bound,
)

def calculate_delG_transport(self):
"""Gibbs energy of transport has two components, 1) proton transport across the membrane, This needs to be adjusted for compartment specific pH 2) charge transport, which is affected by membrane potential.

Expand All @@ -270,7 +325,7 @@ def calculate_delG_transport(self):
Gibbs energy of transport
"""
if len(self.compartments) != 2:
return 0
return 0, 0
else:
charge_dict, proton_dict = self.calculate_transport_charge()
comps = list(charge_dict.keys())
Expand All @@ -291,7 +346,7 @@ def calculate_delG_transport(self):
)
)

return electro_static_delG + proton_potential_adjustment
return electro_static_delG, proton_potential_adjustment

def cal_stoichiometric_matrix(self):
"""stoichiometric vector of forward reaction. Vector has length of number of metabolites in the model.
Expand Down
Loading
Loading