0001 function [A, b_L, b_U, model] = defineLinearConstraints(model, method)
0002
0003 fOffset = 1e-4;
0004
0005 if nargin < 2
0006 method = 1;
0007 end
0008
0009 b_L = model.lb;
0010 b_U = model.ub;
0011 A = model.N;
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064 marked = zeros(size(model.lb));
0065 n = length(model.lb);
0066 dirs = zeros(size(model.lb));
0067 for i = 1:n
0068 if model.lb(i) == 0 && model.ub(i) > fOffset
0069 dirs(i) = 1;
0070 elseif model.lb(i) < -fOffset && model.ub(i) == 0
0071 dirs(i) = -1;
0072 end
0073 end
0074
0075
0076
0077 for i = 1:n
0078 if mod(i,10) == 0
0079 i;
0080 end
0081
0082 LPproblem.A = model.S;
0083
0084 LPproblem.b = zeros(length(model.mets),1);
0085 LPproblem.csense = ['E'*ones(length(model.mets), 1)];
0086 LPproblem.csense = char(LPproblem.csense);
0087 LPproblem.ub = model.ub;
0088 LPproblem.lb = model.lb;
0089 LPproblem.lb(marked > 0) = fOffset;
0090 LPproblem.ub(marked < 0) = -fOffset;
0091
0092 LPproblem.osense = 1;
0093
0094
0095 LPproblem.c = -dirs;
0096 LPproblem.c(marked == 0 ) = 0;
0097 soln = solveCobraLP(LPproblem);
0098 vm = soln.full;
0099
0100 if soln.stat ~= 1
0101 soln
0102 i
0103 pause;
0104 end
0105
0106 [model.lb model.ub dirs vm];
0107
0108 indexes = (vm.*dirs > fOffset) & marked == 0 & dirs ~= 0;
0109 markednew = marked;
0110 markednew(indexes) = dirs(indexes);
0111
0112 [model.lb model.ub dirs vm markednew];
0113
0114 if markednew == marked
0115 break;
0116 else
0117 marked = markednew;
0118 end
0119
0120
0121
0122
0123
0124
0125
0126 end
0127
0128 if method == 1
0129 if isfield(model, 'N')
0130 A = model.N;
0131 else
0132 A = null(model.S);
0133 end
0134 elseif method ==2
0135 A = model.S;
0136 end
0137 b_L(marked > 0) = fOffset;
0138 b_U(marked < 0) = -fOffset;
0139 [b_L b_U];
0140 x = sum(abs(A),2) < 1e-9;
0141 b_L = b_L(~x);
0142 b_U = b_U(~x);
0143 A = A(~x, :);
0144 return