Skip to content

Commit

Permalink
New split variable fix + cleaned up essential_reactions and essential…
Browse files Browse the repository at this point in the history 
…_genes
  • Loading branch information
@phantomas1234
phantomas1234 committed Feb 26, 2015
1 parent 970afdc commit c00cbbc
Show file tree
Hide file tree
Showing 3 changed files with 125 additions and 55 deletions.
37 changes: 22 additions & 15 deletions cameo/core/reaction.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -119,14 +119,15 @@ def reverse_variable(self):

@property
def lower_bound(self):
model = self.model
if model is not None:
if model.reversible_encoding == 'split' and self.reversibility:
return -1 * self.reverse_variable.ub
else:
return self.variable.lb
else:
return self._lower_bound
return self._lower_bound
# model = self.model
# if model is not None:
# if model.reversible_encoding == 'split' and self.reversibility:
# return -1 * self.reverse_variable.ub
# else:
# return self.variable.lb
# else:
# return self._lower_bound

@lower_bound.setter
def lower_bound(self, value):
Expand All @@ -146,7 +147,7 @@ def lower_bound(self, value):
variable = self.variable
reverse_variable = self.reverse_variable

if model.reversible_encoding == 'split' and value < 0 and self._upper_bound >= 0:
if model.reversible_encoding == 'split' and value < 0 and self._upper_bound > 0:
if self._lower_bound > 0:
variable.lb = 0
try:
Expand All @@ -159,16 +160,18 @@ def lower_bound(self, value):
variable.lb = value
except ValueError:
variable.ub = value
self._upper_bound = value
variable.lb = value

self._lower_bound = value

@property
def upper_bound(self):
if self.model is not None:
return self.variable.ub
else:
return self._upper_bound
return self._upper_bound
# if self.model is not None:
# return self.variable.ub
# else:
# return self._upper_bound

@upper_bound.setter
def upper_bound(self, value):
Expand All @@ -181,22 +184,23 @@ def upper_bound(self, value):
reverse_variable.lb, reverse_variable.ub = 0, 0

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

if model.reversible_encoding == 'split' and value > 0 and self._lower_bound < 0:
variable.ub = value
if self._upper_bound < 0:
if self._upper_bound <= 0:
reverse_variable.ub = -1 * variable.lb
variable.lb = 0
else:
try:
variable.ub = value
except ValueError:
variable.lb = value
self._lower_bound = value
variable.ub = value

self._upper_bound = value
Expand Down Expand Up @@ -260,8 +264,11 @@ def add_metabolites(self, metabolites, **kwargs):

def knock_out(self, time_machine=None):
def _(reaction, lb, ub):
# print 'reaction, lb, ub', reaction, lb, ub
# print 'reaction1, reaction.lower_bound, reaction.upper_bound', reaction, reaction.lower_bound, reaction.upper_bound
reaction.upper_bound = ub
reaction.lower_bound = lb
# print 'reaction2, reaction.lower_bound, reaction.upper_bound', reaction, reaction.lower_bound, reaction.upper_bound
if time_machine is not None:
time_machine(do=super(Reaction, self).knock_out, undo=partial(_, self, self.lower_bound, self.upper_bound))
else:
Expand Down
82 changes: 44 additions & 38 deletions cameo/core/solver_based_model.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -421,67 +421,73 @@ def __dir__(self):
return fields

def essential_reactions(self, threshold=1e-6):
"""Return a list of essential reactions.
Parameters
----------
threshold : float (default 1e-6)
Minimal objective flux to be considered viable.
Returns
-------
list
List of essential reactions
"""
essential = []
time_machine = TimeMachine()
try:
solution = self.solve()
except SolveError as e:
print 'Cannot determine essential reactions for un-optimal model.'
raise e
for reaction_id, flux in solution.x_dict.iteritems():
for reaction_id, flux in solution.fluxes.iteritems():
if abs(flux) > 0:
reaction = self.reactions.get_by_id(reaction_id)
# time_machine(do=partial(setattr, reaction, 'lower_bound', 0),
# undo=partial(setattr, reaction, 'lower_bound', reaction.lower_bound))
# time_machine(do=partial(setattr, reaction, 'upper_bound', 0),
# undo=partial(setattr, reaction, 'upper_bound', reaction.upper_bound))
reaction.knock_out(time_machine=time_machine)
try:
sol = self.solve()
except (Infeasible, UndefinedSolution):
print 'except (Infeasible, UndefinedSolution):',reaction_id
essential.append(reaction)
else:
if reaction.id == 'PGK':
print 'PGK', sol.f, sol.f < threshold
if sol.f < threshold:
with TimeMachine() as tm:
reaction.knock_out(time_machine=tm)
try:
sol = self.solve()
except (Infeasible, UndefinedSolution):
essential.append(reaction)
if reaction.id == 'PGK':
print 'PGK asdf', sol.f
print essential
finally:
time_machine.reset()
else:
if sol.f < threshold:
essential.append(reaction)
return essential

def essential_genes(self, threshold=1e-6):
"""Return a list of essential genes.
Parameters
----------
threshold : float (default 1e-6)
Minimal objective flux to be considered viable.
Returns
-------
list
List of essential genes
"""
essential = []
time_machine = TimeMachine()
try:
solution = self.solve()
except SolveError as e:
print 'Cannot determine essential genes for unoptimal model.'
print 'Cannot determine essential genes for un-optimal model.'
raise e
genes_to_check = set()
for reaction_id, flux in solution.x_dict.iteritems():
for reaction_id, flux in solution.fluxes.iteritems():
if abs(flux) > 0:
genes_to_check.update(self.reactions.get_by_id(reaction_id).genes)
for gene in genes_to_check:
reactions = _cobrapy.manipulation.delete.find_gene_knockout_reactions(self, [gene])
for reaction in reactions:
time_machine(do=partial(setattr, reaction, 'lower_bound', 0),
undo=partial(setattr, reaction, 'lower_bound', reaction.lower_bound))
time_machine(do=partial(setattr, reaction, 'upper_bound', 0),
undo=partial(setattr, reaction, 'upper_bound', reaction.upper_bound))
try:
sol = self.solve()
except (Infeasible, UndefinedSolution):
essential.append(gene)
else:
if sol.f < threshold:
with TimeMachine() as tm:
for reaction in reactions:
reaction.knock_out(time_machine=tm)
try:
sol = self.solve()
except (Infeasible, UndefinedSolution):
essential.append(gene)
time_machine.reset()
finally:
time_machine.reset()
else:
if sol.f < threshold:
essential.append(gene)
return essential

def medium(self):
Expand Down
61 changes: 59 additions & 2 deletions tests/test_solver_based_model.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -239,13 +239,15 @@ def test_make_irreversible_irreversible_to_the_other_side(self):
self.assertEqual(pfk_reaction.variable.ub, 999999.)
self.assertEqual(pfk_reaction.reverse_variable, None)
pfk_reaction.upper_bound = -100.
self.assertEqual(pfk_reaction.reverse_variable, None)
pfk_reaction.lower_bound = -999999.
self.assertEqual(pfk_reaction.lower_bound, -999999.)
self.assertEqual(pfk_reaction.upper_bound, -100.)
self.assertEqual(pfk_reaction.variable.lb, -999999.)
self.assertEqual(pfk_reaction.variable.ub, -100)
self.assertEqual(pfk_reaction.reverse_variable.lb, 0.)
self.assertEqual(pfk_reaction.reverse_variable.ub, 0.)
self.assertEqual(pfk_reaction.reverse_variable, None)
# self.assertEqual(pfk_reaction.reverse_variable.lb, 0.)
# self.assertEqual(pfk_reaction.reverse_variable.ub, 0.)

def test_make_lhs_irreversible_reversible(self):
model = self.model
Expand Down Expand Up @@ -302,6 +304,61 @@ def test_knockout(self):
self.assertEqual(reaction.lower_bound, original_bounds[reaction.id][0])
self.assertEqual(reaction.upper_bound, original_bounds[reaction.id][1])

def test_weird_left_to_right_reaction_issue(self):

model = Model("Toy Model")

m1 = Metabolite("M1")
d1 = Reaction("ex1")
d1.add_metabolites({m1: -1})
d1.upper_bound = 0
d1.lower_bound = -1000
# print d1.reaction, d1.lower_bound, d1.upper_bound
model.add_reactions([d1])
self.assertFalse(d1.reversibility)
self.assertEqual(d1.lower_bound, -1000)
self.assertEqual(d1._lower_bound, -1000)
self.assertEqual(d1.upper_bound, 0)
self.assertEqual(d1._upper_bound, 0)
with TimeMachine() as tm:
d1.knock_out(time_machine=tm)
self.assertEqual(d1.lower_bound, 0)
self.assertEqual(d1._lower_bound, 0)
self.assertEqual(d1.upper_bound, 0)
self.assertEqual(d1._upper_bound, 0)
self.assertEqual(d1.lower_bound, -1000)
self.assertEqual(d1._lower_bound, -1000)
self.assertEqual(d1.upper_bound, 0)
self.assertEqual(d1._upper_bound, 0)

def test_one_left_to_right_reaction_set_positive_ub(self):

model = Model("Toy Model")

m1 = Metabolite("M1")
d1 = Reaction("ex1")
d1.add_metabolites({m1: -1})
d1.upper_bound = 0
d1.lower_bound = -1000
model.add_reactions([d1])
self.assertEqual(d1.reverse_variable, None)
self.assertEqual(d1._lower_bound, -1000)
self.assertEqual(d1.lower_bound, -1000)
self.assertEqual(d1._upper_bound, 0)
self.assertEqual(d1.upper_bound, 0)
self.assertEqual(d1.variable.lb, -1000)
self.assertEqual(d1.variable.ub, 0)
d1.upper_bound = .1
self.assertTrue(d1.reverse_variable is not None)
self.assertEqual(d1.variable.lb, 0)
self.assertEqual(d1.variable.ub, .1)
self.assertEqual(d1.reverse_variable.lb, 0)
self.assertEqual(d1.reverse_variable.ub, 1000)
self.assertEqual(d1._lower_bound, -1000)
self.assertEqual(d1.upper_bound, .1)
self.assertEqual(d1._lower_bound, -1000)
self.assertEqual(d1.upper_bound, .1)

def test_iMM904_4HGLSDm_problem(self):
model = load_model(os.path.join(TESTDIR, 'data/iMM904.xml'))
# set upper bound before lower bound after knockout
Expand Down

0 comments on commit c00cbbc

Please sign in to comment.