/
enumerateOptimalSolutions.m
92 lines (82 loc) · 3.43 KB
/
enumerateOptimalSolutions.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
function [solution] = enumerateOptimalSolutions(model)
% Returns a set of optimal flux distributions
% spanning the optimal set
%
% USAGE:
%
% [solution] = enumerateOptimalSolution(model)
%
% INPUT:
% model: COBRA model structure
%
% OUTPUT:
% solution: solution structure
%
% * fluxes - Flux distribution for each iteration
% * nonzero - Boolean matrix denoting which fluxes are nonzero for each iteration
%
% .. Authors:
% - Jan Schellenberger, August 2008 - Based on code by Jennie Reed
% - Reed, J.L. and Palsson, B.O., "Genome-scale in silico models of ''E. coli'' have multiple equivalent phenotypic states: assessment of correlated reaction subsets that comprise network states" , Genome Research, 14:1797-1805(2004).
[nMets,nRxns] = size(model.S);
solution.fluxes = zeros(nRxns,0);
solution.nonzero = zeros(nRxns,0);
%sol = optimizeCbModel(model);
LPproblem = buildLPproblemFromModel(model);
NZ = zeros(nRxns,0);
sol = solveCobraLP(LPproblem);
tol = getCobraSolverParams('LP','feasTol');
NZ(:,end+1) = abs(sol.full(1:nRxns))>tol;
PrevNZ = NZ(:,end);
PrevFW = sol.full(1:nRxns) > tol;
solution.fluxes(:,end+1) = sol.full(1:nRxns);
maxObjective = sol.obj;
MILPproblem = LPproblem;
% Add a Line to keep the objective value
MILPproblem.A = [MILPproblem.A;LPproblem.c'];
MILPproblem.b = [MILPproblem.b;maxObjective];
if MILPproblem.osense == 1
sense = 'L';
else
sense = 'G';
end
MILPproblem.csense = [MILPproblem.csense;sense];
% add the indicator variables to the LP:
[nCtrs,nVars] = size(MILPproblem.A);
% v_+ / v_-
MILPproblem.A = [MILPproblem.A,sparse(nCtrs,nRxns*2);...
speye(nRxns,nVars), -diag(model.ub), sparse(nRxns,nRxns);... %v - v_+*ub <= 0
speye(nRxns,nVars+nRxns), -diag(model.lb);....%v - v_-*lb >= 0
sparse(nRxns,nVars),speye(nRxns,nRxns), speye(nRxns,nRxns)];%;...% v_i+ + v_i- <= 1
%sparse(1,nVars),~(PrevFW & PrevNZ)',~(~PrevFW & PrevNZ)']; %At least one reaction that was not yet active has to be active now!
MILPproblem.vartype = [repmat('C',nVars,1);repmat('B',2*nRxns,1)];
MILPproblem.lb = [MILPproblem.lb;zeros(2*nRxns,1)];
MILPproblem.ub = [MILPproblem.ub;ones(2*nRxns,1)];
MILPproblem.csense = [MILPproblem.csense;...
repmat('L',nRxns,1);...
repmat('G',nRxns,1);...
repmat('L',nRxns,1)];
MILPproblem.b = [MILPproblem.b;...
zeros(2*nRxns,1);...
ones(nRxns,1)];
MILPproblem.c = [LPproblem.c; zeros(2*nRxns,1)];
NonZeroConstPos = size(MILPproblem.A,1);
while 1
% modify the "At least one new reaction" constraint
%MILPproblem.A(NonZeroConstPos,:) = [sparse(1,nVars),~PrevNZ',~PrevNZ'];
% add the "this solution is no longer allowed constraint
MILPproblem.A = [MILPproblem.A;...
sparse(1,nVars), (PrevFW & PrevNZ)', (~PrevFW & PrevNZ)'];
MILPproblem.csense = [MILPproblem.csense; 'L'];
MILPproblem.b = [MILPproblem.b; sum(PrevNZ)-1];
solMILP = solveCobraMILP(MILPproblem);
if solMILP.stat ~= 1
% No more solutions can be found
solution.nonzero = NZ;
return;
end
NZ(:,end+1) = abs(solMILP.full(1:nRxns))>tol;
PrevNZ = NZ(:,end);
PrevFW = solMILP.full(1:nRxns) > tol;
solution.fluxes(:,end+1) = solMILP.full(1:nRxns);
end