Skip to content

Commit

Permalink
efm in progress
Browse files Browse the repository at this point in the history
  • Loading branch information
@the-code-magician
the-code-magician committed Sep 18, 2014
1 parent e01bc6f commit 886e814
Show file tree
Hide file tree
Showing 3 changed files with 59 additions and 51 deletions.
88 changes: 48 additions & 40 deletions cameo/flux_analysis/simulation.py
View file
Expand Up @@ -19,6 +19,10 @@
from sympy import Add
from sympy import Mul

import logging

logger = logging.getLogger("cameo")

add = Add._from_args
mul = Mul._from_args
NegativeOne = sympy.singleton.S.NegativeOne
Expand Down Expand Up @@ -47,6 +51,7 @@ def fba(model, objective=None, *args, **kwargs):
try:
solution = model.solve()
tm.reset()
logger.debug("FBA z=%f" % solution.f)
return solution
except SolveError as e:
raise e
Expand All @@ -71,37 +76,35 @@ def pfba(model, objective=None, *args, **kwargs):
original_objective = model.objective.expression
tm(do=partial(setattr, model, 'reversible_encoding', 'split'),
undo=partial(setattr, model, 'reversible_encoding', model.reversible_encoding))
try:
if objective is not None:
tm(do=partial(setattr, model, 'objective', objective),
undo=partial(setattr, model, 'objective', original_objective))
try:
obj_val = model.solve().f
except SolveError as e:
print "pfba could not determine maximum objective value for\n%s." % model.objective
raise e
if model.objective.direction == 'max':
fix_obj_constraint = model.solver.interface.Constraint(model.objective.expression, lb=obj_val)
else:
fix_obj_constraint = model.solver.interface.Constraint(model.objective.expression, ub=obj_val)
tm(do=partial(model.solver._add_constraint, fix_obj_constraint),
undo=partial(model.solver._remove_constraint, fix_obj_constraint))

pfba_obj = model.solver.interface.Objective(add(
[mul((sympy.singleton.S.One, variable)) for variable in model.solver.variables.values()]),
direction='min', sloppy=True)
# tic = time.time()
tm(do=partial(setattr, model, 'objective', pfba_obj),

if objective is not None:
tm(do=partial(setattr, model, 'objective', objective),
undo=partial(setattr, model, 'objective', original_objective))
# print "obj: ", time.time() - tic
try:
solution = model.solve()
tm.reset()
return solution
except SolveError as e:
print "pfba could not determine an optimal solution for objective %s" % model.objective
raise e
except Exception as e:
try:
obj_val = model.solve().f
except SolveError as e:
print "pfba could not determine maximum objective value for\n%s." % model.objective
raise e
if model.objective.direction == 'max':
fix_obj_constraint = model.solver.interface.Constraint(model.objective.expression, lb=obj_val)
else:
fix_obj_constraint = model.solver.interface.Constraint(model.objective.expression, ub=obj_val)
tm(do=partial(model.solver._add_constraint, fix_obj_constraint),
undo=partial(model.solver._remove_constraint, fix_obj_constraint))

pfba_obj = model.solver.interface.Objective(add(
[mul((sympy.singleton.S.One, variable)) for variable in model.solver.variables.values()]),
direction='min', sloppy=True)
# tic = time.time()
tm(do=partial(setattr, model, 'objective', pfba_obj),
undo=partial(setattr, model, 'objective', original_objective))
# print "obj: ", time.time() - tic
try:
solution = model.solve()
logger.debug("PFBA z=%f" % solution.f)
return solution
except SolveError as e:
print "pfba could not determine an optimal solution for objective %s" % model.objective
raise e


Expand Down Expand Up @@ -170,7 +173,8 @@ def lmoma(model, reference=None, cache={}, volatile=True, *args, **kwargs):
else:
constraint_a = model.solver.interface.Constraint(add([pos_var, mul([NegativeOne, reaction_variable])]),
lb=-flux_value,
sloppy=True)
sloppy=True,
name=constraint_a_id)
if not volatile:
cache['constraints'][constraint_a_id] = constraint_a

Expand All @@ -181,7 +185,8 @@ def lmoma(model, reference=None, cache={}, volatile=True, *args, **kwargs):
else:
constraint_b = model.solver.interface.Constraint(add([neg_var, reaction.variable]),
lb=flux_value,
sloppy=True)
sloppy=True,
name=constraint_b_id)
if not volatile:
cache['constraints'][constraint_b_id] = constraint_b

Expand All @@ -197,6 +202,7 @@ def lmoma(model, reference=None, cache={}, volatile=True, *args, **kwargs):

try:
solution = model.solve()
logger.debug("LMOMA z=%f" % solution.f)
return solution
except SolveError as e:
#print "lmoma could not determine an optimal solution for objective %s" % model.objective
Expand Down Expand Up @@ -258,12 +264,10 @@ def room(model, reference=None, cache={}, volatile=True, delta=0.03, epsilon=0.0
constraint_a.ub = w_u
else:

# vi - yi(vmaxi + w_ui) >= w_ui
expression = add([
reaction.variable,
mul([RealNumber(-(reaction.upper_bound + w_u)), var])])
# vi - yi(vmaxi + w_ui) <= w_ui
expression = add([reaction.variable, mul([RealNumber(-(reaction.upper_bound + w_u)), var])])

constraint_a = (model.solver.interface.Constraint(expression, ub=w_u, sloppy=True))
constraint_a = model.solver.interface.Constraint(expression, ub=w_u, name=constraint_a_id, sloppy=True)
if not volatile:
cache['constraints'][constraint_a_id] = constraint_a

Expand All @@ -276,12 +280,12 @@ def room(model, reference=None, cache={}, volatile=True, delta=0.03, epsilon=0.0
constraint_b._set_coefficients_low_level({var: -reaction.lower_bound + w_l})
constraint_b.lb = w_l
else:
# vi - yi(vmini - w_li) <= w_li
# vi - yi(vmini - w_li) >= w_li
expression = add([
reaction.variable,
mul([RealNumber(-(reaction.lower_bound - w_l)), var])])

constraint_b = (model.solver.interface.Constraint(expression, lb=w_l, sloppy=True))
constraint_b = model.solver.interface.Constraint(expression, lb=w_l, name=constraint_b_id, sloppy=True)
if not volatile:
cache['constraints'][constraint_b_id] = constraint_b

Expand All @@ -297,7 +301,9 @@ def room(model, reference=None, cache={}, volatile=True, delta=0.03, epsilon=0.0
cache['first_run'] = False

try:
return model.solve()
solution = model.solve()
logger.debug("ROOM z=%f" % solution.f)
return solution
except SolveError as e:
print "room could not determine an optimal solution for objective %s" % model.objective
raise e
Expand Down Expand Up @@ -395,7 +401,9 @@ def reset_model(model, cache):
from cobra.io import read_sbml_model
from cobra.flux_analysis.parsimonious import optimize_minimal_flux
from cameo import load_model
import logging

logging.getLogger("cameo").setLevel(logging.DEBUG)
# sbml_path = '../../tests/data/EcoliCore.xml'
sbml_path = '../../tests/data/iJO1366.xml'

Expand Down
3 changes: 1 addition & 2 deletions cameo/io.py
View file
Expand Up @@ -64,5 +64,4 @@ def load_model(path_or_handle, solver_interface=optlang.glpk_interface):
else:
if model.interface is not solver_interface and solver_interface is not None:
model.solver = solver_interface
return model

return model
19 changes: 10 additions & 9 deletions cameo/structural_analysis/__init__.py
View file
Expand Up @@ -97,15 +97,15 @@ class EFMModel(object):
.. [2] Axel von Kamp, Steffen Klamt "Enumeration of Smallest Intervention Strategies in Genome-Scale
Metabolic Networks" PLOS Computational Biology, vol 10, issue 01, pp. e1003378, 2014.
"""
def __init__(self, model, M=100000, matrix=None, include_exchanges=False):
def __init__(self, model, M=100000, matrix=None, include_exchanges=False, flux_type="integer"):
self.z_map = dict()
self.t_map = dict()
self.M = M
if matrix is None:
matrix = model.S
self.matrix = matrix
self.model = model.solver.interface.Model()
self._populate_model(model, include_exchanges)
self._populate_model(model, include_exchanges, flux_type)

def add(self, obj):
self.model.add(obj)
Expand All @@ -127,14 +127,14 @@ def t_vars(self):
def t_var(self, reaction_id):
return self.t_map[reaction_id]

def _populate_model(self, model, include_exchanges):
def _populate_model(self, model, include_exchanges, flux_type):
constraints = list()
exchanges = model.exchanges
for reaction in model.reactions:
if not reaction in exchanges or include_exchanges:
z, t = self.add_reaction(reaction.id, model, constraints)
z, t = self.add_reaction(reaction.id, model, constraints, flux_type)
if reaction.reversibility:
z_rev, t_rev = self.add_reaction(reaction._get_reverse_id(), model, constraints)
z_rev, t_rev = self.add_reaction(reaction._get_reverse_id(), model, constraints, flux_type)
exp = add([z, z_rev])
rev_constraint = model.solver.interface.Constraint(exp, ub=1, name="rev_%s" % reaction.id)
constraints.append(rev_constraint)
Expand Down Expand Up @@ -168,10 +168,10 @@ def _populate_model(self, model, include_exchanges):
self.model.add(constraints)
self.model.objective = objective

def add_reaction(self, r_id, model, constraints):
def add_reaction(self, r_id, model, constraints, flux_type):
z = model.solver.interface.Variable("z_%s" % r_id, type='binary')
self.z_map[r_id] = z
t = model.solver.interface.Variable("t_%s" % r_id, lb=0, type='integer')
t = model.solver.interface.Variable("t_%s" % r_id, lb=0, ub=self.M, type=flux_type)
self.t_map[r_id] = t

expression_1 = add([mul([RealNumber(-self.M), z]), t])
Expand Down Expand Up @@ -216,6 +216,7 @@ def shortest_elementary_flux_modes(model=None, k=5, M=100000, efm_model=None, ma
genome-scale metabolic networks" Bioinformatics, vol 25, issue 23, pp. 3158-3165, 2009.
"""
if matrix is None:
logger.debug("Computing integer matrix...")
matrix = model.integerS

if efm_model is None:
Expand Down Expand Up @@ -247,7 +248,7 @@ def shortest_elementary_flux_modes(model=None, k=5, M=100000, efm_model=None, ma
efm_model.add(iteration_constraint)

#interation n
logger.debug("Iteration %i:" % i+2)
logger.debug("Iteration %i:" % (i+2))
status = efm_model.optimize()
logger.debug("Iteration %i: %s" % (i+2, status))
[logger.debug("%s: %f" % (z, z.primal)) for z in efm_model.z_vars if z.primal > 0]
Expand All @@ -270,7 +271,7 @@ def fixed_size_elementary_modes(model=None, c=1, efm_model=None, size=20):
"""

if efm_model is None:
efm_model = EFMModel(model, M=c, include_exchanges=True, matrix=model.S)
efm_model = EFMModel(model, M=c, include_exchanges=True, matrix=model.S, flux_type="continuous")

with TimeMachine() as tm:
expression = sum(efm_model.z_vars)
Expand Down

0 comments on commit 886e814

Please sign in to comment.