Skip to content

Commit

Permalink
Merge pull request #27 from biosustain/reverse_variables_for_everyone
Browse files Browse the repository at this point in the history 
All reactions (not just reversible ones) have now forward_variable and r...
  • Loading branch information
@phantomas1234
phantomas1234 committed Apr 10, 2015
2 parents a7195a7 + 285b2df commit d1d64c6
Show file tree
Hide file tree
Showing 7 changed files with 193 additions and 184 deletions.
176 changes: 113 additions & 63 deletions cameo/core/reaction.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,10 +15,11 @@

from __future__ import absolute_import, print_function

import warnings

from functools import partial
import hashlib
import cobra as _cobrapy
import sympy

import cameo
from cameo import flux_analysis
Expand Down Expand Up @@ -88,22 +89,46 @@ def __str__(self):
return ''.join((self.id, ": ", self.build_reaction_string()))

@property
def variable(self):
def reversibility(self):
return self._lower_bound < 0 and self._upper_bound > 0

def _get_reverse_id(self):
"""Generate the id of reverse_variable from the reaction's id."""
return '_'.join((self.id, 'reverse', hashlib.md5(self.id.encode('utf-8')).hexdigest()[0:5]))

def _get_forward_id(self):
"""Generate the id of forward_variable from the reaction's id."""
return self.id
# return '_'.join((self.id, 'forward', hashlib.md5(self.id.encode('utf-8')).hexdigest()[0:5]))

@property
def flux_expression(self):
"""An optlang variable representing the forward flux (if associated with model), otherwise None.
Representing the net flux if model.reversible_encoding == 'unsplit'"""
model = self.model
if model is not None:
return model.solver.variables[self.id]
return 1.*self.forward_variable - 1.*self.reverse_variable
else:
return None

@property
def reversibility(self):
return self._lower_bound < 0 and self._upper_bound > 0
def variable(self):
warnings.warn('reaction.variable is deprecated. Please use reaction.forward_variable.', DeprecationWarning)
return self.forward_variable

def _get_reverse_id(self):
"""Generate the id of revers_variable from the reaction's id."""
return '_'.join((self.id, 'reverse', hashlib.md5(self.id.encode('utf-8')).hexdigest()[0:5]))
@property
def forward_variable(self):
"""An optlang variable representing the forward flux (if associated with model), otherwise None.
Representing the net flux if model.reversible_encoding == 'unsplit'"""
model = self.model
if model is not None:
aux_id = self._get_forward_id()
try:
return model.solver.variables[aux_id]
except KeyError:
return None
else:
return None

@property
def reverse_variable(self):
Expand All @@ -128,33 +153,48 @@ def lower_bound(self, value):

if model is not None:

if value >= 0 and self._lower_bound < 0 and self._upper_bound > 0:
reverse_variable = self.reverse_variable
reverse_variable.lb, reverse_variable.ub = 0, 0
elif value < 0 and self._lower_bound >= 0 and self.reverse_variable is None: # self._lower_bound >= 0 implies self._upper_bound >= 0
reverse_variable = model.solver._add_variable(
model.solver.interface.Variable(self._get_reverse_id(), lb=0, ub=0))
for met, coeff in six.iteritems(self._metabolites):
model.solver.constraints[met.id] += sympy.Mul._from_args((-1 * sympy.RealNumber(coeff), reverse_variable))

variable = self.variable
reverse_variable = self.reverse_variable

if model.reversible_encoding == 'split' and value < 0 and self._upper_bound > 0:
if self._lower_bound >= 0:
variable.lb = 0
try:
reverse_variable.ub = -1 * value
except ValueError:
reverse_variable.lb = -1 * value
reverse_variable.ub = -1 * value
else:
try:
variable.lb = value
except ValueError:
variable.ub = value
forward_variable, reverse_variable = self.forward_variable, self.reverse_variable
if self._lower_bound < 0 and self._upper_bound > 0: # reversible
if value < 0:
reverse_variable.ub = -1*value
elif value >= 0:
reverse_variable.ub = 0
try:
forward_variable.lb = value
except ValueError:
forward_variable.ub = value
self._upper_bound = value
forward_variable.lb = value
elif self._lower_bound == 0 and self._upper_bound == 0: # knockout
if value < 0:
reverse_variable.ub = -1*value
elif value >= 0:
forward_variable.ub = value
forward_variable.lb = value
elif self._lower_bound >= 0: # forward irreversible
if value < 0:
reverse_variable.ub = -1*value
forward_variable.lb = 0
else:
try:
forward_variable.lb = value
except ValueError:
forward_variable.ub = value
self._upper_bound = value
forward_variable.lb = value

elif self._upper_bound <= 0: # reverse irreversible
if value > 0:
reverse_variable.ub = 0
reverse_variable.ub = 0
forward_variable.ub = value
self._upper_bound = value
variable.lb = value
forward_variable.lb = value
else:
reverse_variable.ub = -1*value
else:
print({'value': value, 'self._lower_bound': self._lower_bound, 'self._upper_bound': self._upper_bound})
raise ValueError('lower_bound issue')

self._lower_bound = value

Expand All @@ -166,38 +206,48 @@ def upper_bound(self):
def upper_bound(self, value):
model = self.model
if model is not None:
# Remove auxiliary variable if not needed anymore
reverse_variable = self.reverse_variable
variable = self.variable
if value <= 0 and self._upper_bound > 0 and self._lower_bound < 0:
reverse_variable.lb, reverse_variable.ub = 0, 0
try:
variable.ub = value
variable.lb = self._lower_bound
except ValueError:
variable.lb = value

forward_variable, reverse_variable = self.forward_variable, self.reverse_variable
if self._lower_bound < 0 and self._upper_bound > 0: # reversible
if value > 0:
forward_variable.ub = value
elif value <= 0:
forward_variable.ub = 0
try:
reverse_variable.lb = -1*value
except ValueError:
reverse_variable.ub = -1*value
self._lower_bound = value
reverse_variable.lb = -1*value
elif self._lower_bound == 0 and self._upper_bound == 0: # knockout
if value > 0:
forward_variable.ub = value
elif value <= 0:
reverse_variable.ub = -1*value
elif self._lower_bound >= 0: # forward irreversible
if value > 0:
forward_variable.ub = value
else:
forward_variable.lb = 0
forward_variable.ub = 0
reverse_variable.ub = -1*value
self._lower_bound = value
variable.ub = value

# Add auxiliary variable if needed
elif value > 0 and self._upper_bound <= 0: # self._upper_bound < 0 implies self._lower_bound < 0
if self.reverse_variable is None:
reverse_variable = model.solver._add_variable(
model.solver.interface.Variable(self._get_reverse_id(), lb=0, ub=0))
for met, coeff in six.iteritems(self._metabolites):
model.solver.constraints[met.id] += sympy.Mul._from_args((-1 * sympy.RealNumber(coeff), reverse_variable))
reverse_variable.lb = -1*value

elif self._upper_bound <= 0: # reverse irreversible
if value < 0:
try:
reverse_variable.lb = -1*value
except ValueError:
reverse_variable.ub = -1*value
self._lower_bound = value
reverse_variable.lb = -1*value
else:
reverse_variable = self.reverse_variable
variable.ub = value
reverse_variable.ub = -1 * variable.lb
variable.lb = 0
forward_variable.ub = value
reverse_variable.lb = 0
else:
try:
variable.ub = value
except ValueError:
variable.lb = value
self._lower_bound = value
variable.ub = value
print({'value': value, 'self._lower_bound': self._lower_bound, 'self._upper_bound': self._upper_bound})
raise ValueError('upper_bound issue')

self._upper_bound = value

Expand Down
67 changes: 22 additions & 45 deletions cameo/core/solver_based_model.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -73,7 +73,6 @@ class SolverBasedModel(_cobrapy.core.Model):

def __init__(self, description=None, solver_interface=optlang, **kwargs):
super(SolverBasedModel, self).__init__(description, **kwargs)
self._reversible_encoding = 'split'
cleaned_reactions = _cobrapy.core.DictList()
for reaction in self.reactions:
if isinstance(reaction, Reaction):
Expand Down Expand Up @@ -147,12 +146,7 @@ def objective(self, value):
if isinstance(value, str):
value = self.reactions.get_by_id(value)
if isinstance(value, Reaction):
if self.reversible_encoding == 'split' and value.reverse_variable is not None:
obj_expression = Add._from_args((Mul._from_args((S.One, value.variable)), Mul._from_args((S.NegativeOne, value.reverse_variable))))
self.solver.objective = self.solver.interface.Objective(obj_expression, sloppy=True)
else:
obj_expression = Mul._from_args((S.One, value.variable))
self.solver.objective = self.solver.interface.Objective(obj_expression, sloppy=True)
self.solver.objective = self.solver.interface.Objective(value.flux_expression, sloppy=True)
elif isinstance(value, self.solver.interface.Objective):
self.solver.objective = value
# TODO: maybe the following should be allowed
Expand Down Expand Up @@ -192,12 +186,13 @@ def _populate_solver_from_scratch(self):
for rxn in self.reactions:
lower_bound, upper_bound = rxn._lower_bound, rxn._upper_bound
if lower_bound < 0 and upper_bound > 0: # i.e. lower_bound < 0 and upper_bound > 0
var = self.solver._add_variable(self.solver.interface.Variable(rxn.id, lb=0, ub=upper_bound))
var = self.solver._add_variable(self.solver.interface.Variable(rxn._get_forward_id(), lb=0, ub=upper_bound))
aux_var = self.solver._add_variable(
self.solver.interface.Variable(rxn._get_reverse_id(), lb=0, ub=-1 * lower_bound))
else:
var = self.solver._add_variable(
self.solver.interface.Variable(rxn.id, lb=lower_bound, ub=upper_bound))
var = self.solver._add_variable(self.solver.interface.Variable(rxn._get_forward_id(), lb=lower_bound, ub=upper_bound))
aux_var = self.solver._add_variable(
self.solver.interface.Variable(rxn._get_reverse_id(), lb=0, ub=0))
if rxn.objective_coefficient != 0.:
objective_terms.append(sympy.Mul._from_args((sympy.RealNumber(rxn.objective_coefficient), var)))
for met, coeff in six.iteritems(rxn._metabolites):
Expand All @@ -220,29 +215,6 @@ def _populate_solver_from_scratch(self):
objective_expression = sympy.Add._from_args(objective_terms)
self.solver.objective = self.solver.interface.Objective(objective_expression, name='obj', direction='max')

@property
def reversible_encoding(self):
return self._reversible_encoding

@reversible_encoding.setter
def reversible_encoding(self, value):
if self._reversible_encoding == value:
pass
else:
if value == 'unsplit':
for reaction in self.reactions:
if reaction.reversibility:
reaction.variable.lb = -1 * reaction.reverse_variable.ub
reaction.reverse_variable.ub = 0
elif value == 'split':
for reaction in self.reactions:
if reaction.reversibility:
reaction.reverse_variable.ub = -1 * reaction.variable.lb
reaction.variable.lb = 0
else:
raise ValueError('%s is not a valid encoding. Try one of %s instead.' % (value, ('unsplit', 'split')))
self._reversible_encoding = value

def add_metabolites(self, metabolite_list):
super(SolverBasedModel, self).add_metabolites(metabolite_list)
for met in metabolite_list:
Expand All @@ -263,24 +235,29 @@ def add_reactions(self, reaction_list):
constr_terms = dict()
metabolites = {}
for reaction in cloned_reaction_list:
if reaction.reversibility and self._reversible_encoding == "split":
reaction_variable = self.solver.interface.Variable(reaction.id, lb=0, ub=reaction._upper_bound)
aux_var = self.solver.interface.Variable(reaction._get_reverse_id(), lb=0, ub=-reaction._lower_bound)
self.solver._add_variable(aux_var)
else:
reaction_variable = self.solver.interface.Variable(reaction.id, lb=reaction._lower_bound,
ub=reaction._upper_bound)
self.solver._add_variable(reaction_variable)
if reaction.reversibility:
forward_variable = self.solver.interface.Variable(reaction._get_forward_id(), lb=0, ub=reaction._upper_bound)
reverse_variable = self.solver.interface.Variable(reaction._get_reverse_id(), lb=0, ub=-1*reaction._lower_bound)
elif 0 == reaction.lower_bound and reaction.upper_bound == 0:
forward_variable = self.solver.interface.Variable(reaction._get_forward_id(), lb=0, ub=0)
reverse_variable = self.solver.interface.Variable(reaction._get_reverse_id(), lb=0, ub=0)
elif reaction.lower_bound >= 0:
forward_variable = self.solver.interface.Variable(reaction.id, lb=reaction._lower_bound, ub=reaction._upper_bound)
reverse_variable = self.solver.interface.Variable(reaction._get_reverse_id(), lb=0, ub=0)
elif reaction.upper_bound <= 0:
forward_variable = self.solver.interface.Variable(reaction.id, lb=0, ub=0)
reverse_variable = self.solver.interface.Variable(reaction._get_reverse_id(), lb=-1*reaction._upper_bound, ub=-1*reaction._lower_bound)
self.solver._add_variable(forward_variable)
self.solver._add_variable(reverse_variable)

for metabolite, coeff in six.iteritems(reaction.metabolites):
if metabolite.id in constr_terms:
constr_terms[metabolite.id].append(
sympy.Mul._from_args([sympy.RealNumber(coeff), reaction_variable]))
sympy.Mul._from_args([sympy.RealNumber(coeff), forward_variable]))
else:
constr_terms[metabolite.id] = [sympy.Mul._from_args([sympy.RealNumber(coeff), reaction_variable])]
constr_terms[metabolite.id] = [sympy.Mul._from_args([sympy.RealNumber(coeff), forward_variable])]
metabolites[metabolite.id] = metabolite
if reaction.reversibility and self._reversible_encoding == "split":
constr_terms[metabolite.id].append(sympy.Mul._from_args([sympy.RealNumber(-1*coeff), aux_var]))
constr_terms[metabolite.id].append(sympy.Mul._from_args([sympy.RealNumber(-1*coeff), reverse_variable]))

for met_id, terms in six.iteritems(constr_terms):
expr = sympy.Add._from_args(terms)
Expand Down
6 changes: 0 additions & 6 deletions cameo/flux_analysis/analysis.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -60,9 +60,6 @@ def flux_variability_analysis(model, reactions=None, fraction_of_optimum=0., rem
if reactions is None:
reactions = model.reactions
with TimeMachine() as tm:
if model.reversible_encoding == 'split':
tm(do=partial(setattr, model, 'reversible_encoding', 'unsplit'),
undo=partial(setattr, model, 'reversible_encoding', 'split'))
if fraction_of_optimum > 0.:
try:
obj_val = model.solve().f
Expand Down Expand Up @@ -117,9 +114,6 @@ def phenotypic_phase_plane(model, variables=[], objective=None, points=20, view=
if view is None:
view = config.default_view
with TimeMachine() as tm:
if model.reversible_encoding == 'split':
tm(do=partial(setattr, model, 'reversible_encoding', 'unsplit'),
undo=partial(setattr, model, 'reversible_encoding', 'split'))
if objective is not None:
tm(do=partial(setattr, model, 'objective', objective),
undo=partial(setattr, model, 'objective', model.objective))
Expand Down
Loading

0 comments on commit d1d64c6

Please sign in to comment.