0001 function [solution] = enumerateOptimalSolutions(model)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021 [m,n] = size(model.S);
0022
0023 solution.fluxes = zeros(n,0);
0024 solution.nonzero = zeros(n,0);
0025
0026 sol = optimizeCbModel(model);
0027 maxObjective = sol.f;
0028
0029
0030 prevNZ = abs(sol.x) > .0001;
0031 NZ = prevNZ;
0032 solution.fluxes = sol.x;
0033 solution.nonzero = prevNZ;
0034
0035 while 1
0036
0037
0038
0039
0040
0041 A = [model.S, zeros(m,2*n)];
0042 b = zeros(m,1);
0043 csense = '';
0044 for i = 1:m
0045 csense(end+1) = 'E';
0046 end
0047
0048 A = [A;
0049 [eye(n,2*n), -diag(model.ub)] ];
0050 b = [b;
0051 zeros(n,1)];
0052 for i = 1:n
0053 csense(end+1) = 'L';
0054 end
0055
0056 A = [A;
0057 eye(n,2*n), -diag(model.lb) ];
0058 b = [b;
0059 zeros(n,1)];
0060 for i = 1:n
0061 csense(end+1) = 'G';
0062 end
0063
0064
0065
0066 A = [A;
0067 zeros(n,n), eye(n,n), eye(n,n) ];
0068 b = [b;
0069 ones(n,1)];
0070 for i = 1:n
0071 csense(end+1) = 'L';
0072 end
0073
0074
0075 A = [A;
0076 zeros(1,n), prevNZ', zeros(1,n) ];
0077 b = [b;
0078 1];
0079 csense(end+1) = 'G';
0080
0081
0082 for i = 1:size(NZ,2)
0083 A = [A;
0084 [zeros(1,n), zeros(1,n) NZ(:,i)']];
0085 b(end+1) = sum(NZ(:,i))-1;
0086 csense(end+1) = 'L';
0087 end
0088
0089
0090 vartype = [];
0091 for i = 1:n
0092 vartype(i) = 'C';
0093 end
0094 for i = 1:2*n
0095 vartype(end+1) = 'B';
0096 end
0097
0098
0099 lb = [model.lb; zeros(2*n,1)];
0100 ub = [model.ub; ones(2*n,1)];
0101
0102 c = [model.c; zeros(2*n,1)];
0103
0104
0105
0106 MILPproblem.A = A;
0107 MILPproblem.b = b;
0108 MILPproblem.c = c;
0109 MILPproblem.csense = csense;
0110 MILPproblem.lb = lb;
0111 MILPproblem.ub = ub;
0112 MILPproblem.osense = -1;
0113 MILPproblem.vartype = vartype;
0114 MILPproblem.x0 = [];
0115
0116
0117
0118
0119 MILPsol = solveCobraMILP(MILPproblem)
0120
0121 NZ(:,end+1) = abs(MILPsol.full(1:n))>.000000001;
0122 PrevNZ = NZ(:,end);
0123
0124
0125 if (abs(MILPsol.full - maxObjective) > .001)
0126 'done';
0127 return;
0128 end
0129 solution.fluxes = [solution.fluxes,MILPsol.full(1:n)];
0130 solution.nonzero = [solution.nonzero, NZ(:,end)];
0131 solution
0132 end
0133
0134 return;