0001 function [MILPproblem] = addLoopLawConstraints(LPproblem, model, rxnIndex)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031 method = 2;
0032 reduce_vars = 1;
0033 combine_vars = 0;
0034
0035 if nargin < 3
0036 if size(LPproblem.A,2) == size(model.S,2);
0037 rxnIndex = 1:size(model.S,2);
0038 elseif size(LPproblem.A,2) > size(model.S,2)
0039 display('warning: extra variables in LPproblem. will assume first n correspond to v')
0040 rxnIndex = 1:size(model.S,2);
0041 else
0042 display('LPproblem must have at least as many variables as model has reactions');
0043 return;
0044 end
0045 elseif length(find(rxnIndex)) ~= size(model.S,2)
0046 display('rxnIndex must contain exactly n entries');
0047 return;
0048 end
0049 if any(rxnIndex > size(LPproblem.A,2))
0050 display('rxnIndex out of bounds');
0051 return;
0052 end
0053
0054 MILPproblem = LPproblem;
0055
0056 S = model.S;
0057 [m,n] = size(LPproblem.A);
0058 nontransport = (sum(S ~= 0) > 1)';
0059
0060 nontransport = (nontransport | (model.lb ==0 & model.ub == 0));
0061
0062
0063 if reduce_vars == 1
0064 active1 = ~(model.lb ==0 & model.ub == 0);
0065 S2 = S(:,active1);
0066
0067 N2 = sparseNull(sparse(S2));
0068 N = zeros(length(active1), size(N2,2));
0069 N(active1,:) = N2;
0070
0071 active = any(abs(N) > 1e-6, 2);
0072
0073
0074 nontransport = nontransport & active;
0075 end
0076
0077 Sn = S(:,nontransport);
0078
0079 Ninternal = sparseNull(sparse(Sn));
0080
0081
0082 linternal = size(Ninternal,2);
0083
0084 nint = length(find(nontransport));
0085 temp = sparse(nint, n);
0086 temp(:, rxnIndex(nontransport)) = speye(nint);
0087
0088
0089 if method == 1
0090 MILPproblem.A = [LPproblem.A, sparse(m,3*nint);
0091 temp, -10000*speye(nint), sparse(nint, 2*nint);
0092 temp, sparse(nint, nint), 10000*speye(nint), sparse(nint, nint);
0093 sparse(nint, n), speye(nint), speye(nint), sparse(nint, nint);
0094 sparse(nint, n), -100*speye(nint), 1*speye(nint), speye(nint);
0095 sparse(nint, n), -1*speye(nint), 100*speye(nint), speye(nint);
0096 sparse(linternal, n+2*nint), Ninternal'];
0097
0098 MILPproblem.b = [LPproblem.b;
0099 zeros(2*nint,1);
0100 ones(nint,1);
0101 zeros(2*nint + linternal,1);];
0102
0103 MILPproblem.c = [LPproblem.c;
0104 zeros(3*nint,1)];
0105
0106 MILPproblem.csense = LPproblem.csense;
0107 for i = 1:nint, MILPproblem.csense(end+1,1) = 'L';end
0108 for i = 1:nint, MILPproblem.csense(end+1,1) = 'G';end
0109 for i = 1:nint, MILPproblem.csense(end+1,1) = 'L';end
0110 for i = 1:nint, MILPproblem.csense(end+1,1) = 'L';end
0111 for i = 1:nint, MILPproblem.csense(end+1,1) = 'G';end
0112 for i = 1:linternal, MILPproblem.csense(end+1,1) = 'E';end
0113
0114 MILPproblem.vartype = [];
0115 if isfield(LPproblem, 'vartype')
0116 MILPproblem.vartype = LPproblem.vartype;
0117 else
0118 for i = 1:n, MILPproblem.vartype(end+1,1) = 'C';end;
0119 end
0120 for i = 1:2*nint, MILPproblem.vartype(end+1,1) = 'B';end;
0121 for i = 1:nint, MILPproblem.vartype(end+1,1) = 'C';end;
0122
0123 if isfield(LPproblem, 'F')
0124 MILPproblem.F = sparse(size(MILPproblem.A,2), size(MILPproblem.A,2));
0125 MILPproblem.F(1:size(LPproblem.F,1), 1:size(LPproblem.F,1)) = LPproblem.F;
0126 end
0127
0128
0129 MILPproblem.lb = [LPproblem.lb;
0130 zeros(nint*2,1);
0131 -1000*ones(nint,1);];
0132 MILPproblem.ub = [LPproblem.ub;
0133 ones(nint*2,1);
0134 1000*ones(nint,1);];
0135
0136 MILPproblem.x0 = [];
0137
0138 elseif method == 2
0139 MILPproblem.A = [LPproblem.A, sparse(m,2*nint);
0140 temp, -10000*speye(nint), sparse(nint, nint);
0141 temp, -10000*speye(nint), sparse(nint, nint);
0142 sparse(nint, n), -101*speye(nint), speye(nint);
0143 sparse(nint, n), -101*speye(nint), speye(nint);
0144 sparse(linternal, n + nint), Ninternal'];
0145
0146 MILPproblem.b = [LPproblem.b;
0147 zeros(nint,1);
0148 -10000*ones(nint, 1);
0149 -ones(nint,1);
0150 -100*ones(nint, 1);
0151 zeros(linternal,1);];
0152
0153 MILPproblem.c = [LPproblem.c;
0154 zeros(2*nint,1)];
0155
0156 MILPproblem.csense = LPproblem.csense;
0157 for i = 1:nint, MILPproblem.csense(end+1,1) = 'L';end
0158 for i = 1:nint, MILPproblem.csense(end+1,1) = 'G';end
0159 for i = 1:nint, MILPproblem.csense(end+1,1) = 'L';end
0160 for i = 1:nint, MILPproblem.csense(end+1,1) = 'G';end
0161 for i = 1:linternal, MILPproblem.csense(end+1,1) = 'E';end
0162
0163 MILPproblem.vartype = '';
0164 if isfield(LPproblem, 'vartype')
0165 MILPproblem.vartype = LPproblem.vartype;
0166 else
0167 for i = 1:n, MILPproblem.vartype(end+1,1) = 'C';end;
0168 end
0169 for i = 1:nint, MILPproblem.vartype(end+1,1) = 'B';end;
0170 for i = 1:nint, MILPproblem.vartype(end+1,1) = 'C';end;
0171
0172 if isfield(LPproblem, 'F')
0173 MILPproblem.F = sparse(size(MILPproblem.A,2), size(MILPproblem.A,2));
0174 MILPproblem.F(1:size(LPproblem.F,1), 1:size(LPproblem.F,1)) = LPproblem.F;
0175 end
0176
0177
0178 MILPproblem.lb = [LPproblem.lb;
0179 zeros(nint,1);
0180 -1000*ones(nint,1);];
0181 MILPproblem.ub = [LPproblem.ub;
0182 ones(nint,1);
0183 1000*ones(nint,1);];
0184
0185 MILPproblem.x0 = [];
0186 elseif method == 3
0187 MILPproblem.A = [LPproblem.A, sparse(m,2*nint);
0188 temp, -10000*speye(nint), sparse(nint, nint);
0189
0190 sparse(nint, n), -101*speye(nint), speye(nint);
0191
0192 sparse(linternal, n + nint), Ninternal'];
0193
0194 MILPproblem.b_L = [LPproblem.b;
0195
0196 -10000*ones(nint, 1);
0197
0198 -100*ones(nint, 1);
0199 zeros(linternal,1);];
0200 MILPproblem.b_U = [LPproblem.b;
0201 zeros(nint,1);
0202
0203 -ones(nint,1);
0204
0205 zeros(linternal,1);];
0206
0207 MILPproblem.b_L(find(LPproblem.csense == 'E')) = LPproblem.b(LPproblem.csense == 'E');
0208 MILPproblem.b_U(find(LPproblem.csense == 'E')) = LPproblem.b(LPproblem.csense == 'E');
0209 MILPproblem.b_L(find(LPproblem.csense == 'G')) = LPproblem.b(LPproblem.csense == 'G');
0210 MILPproblem.b_U(find(LPproblem.csense == 'G')) = inf;
0211 MILPproblem.b_L(find(LPproblem.csense == 'L')) = -inf;
0212 MILPproblem.b_U(find(LPproblem.csense == 'L')) = LPproblem.b(LPproblem.csense == 'L');
0213
0214 MILPproblem.c = [LPproblem.c;
0215 zeros(2*nint,1)];
0216
0217 MILPproblem.csense = [];
0218
0219 MILPproblem.vartype = [];
0220 if isfield(LPproblem, 'vartype')
0221 MILPproblem.vartype = LPproblem.vartype;
0222 else
0223 for i = 1:n, MILPproblem.vartype(end+1,1) = 'C';end;
0224 end
0225 for i = 1:nint, MILPproblem.vartype(end+1,1) = 'B';end;
0226 for i = 1:nint, MILPproblem.vartype(end+1,1) = 'C';end;
0227
0228 if isfield(LPproblem, 'F')
0229 MILPproblem.F = sparse(size(MILPproblem.A,2), size(MILPproblem.A,2));
0230 MILPproblem.F(1:size(LPproblem.F,1), 1:size(LPproblem.F,1)) = LPproblem.F;
0231 end
0232
0233
0234 MILPproblem.lb = [LPproblem.lb;
0235 zeros(nint,1);
0236 -1000*ones(nint,1);];
0237 MILPproblem.ub = [LPproblem.ub;
0238 ones(nint,1);
0239 1000*ones(nint,1);];
0240
0241 MILPproblem.x0 = [];
0242 else
0243 display('method not found')
0244 method
0245 pause;
0246 end
0247
0248 if combine_vars && method == 2
0249
0250
0251 Ns = N(nontransport,:);
0252
0253
0254
0255
0256 Ns2 = Ns;
0257 for i = 1:size(Ns2,1)
0258 m = sqrt(Ns2(i,:)*Ns2(i,:)');
0259 Ns2(i,:) = Ns2(i,:)/m;
0260 end
0261
0262 t = Ns2 * Ns2';
0263
0264
0265
0266
0267
0268
0269
0270
0271
0272 cutoff = .9999999;
0273
0274
0275
0276
0277
0278
0279
0280
0281
0282 t2 = sparse(size(t,1), size(t, 2));
0283 t2(abs(t) > cutoff) = t(abs(t) > cutoff);
0284 t = t2;
0285 checkedbefore = zeros(nint,1);
0286
0287 for i = 2:nint
0288 x = find(t(i,1:i-1)>cutoff);
0289 if ~isempty(x)
0290 checkedbefore(i) = x(1);
0291 end
0292 y = find(t(i,1:i-1)<-cutoff);
0293 if ~isempty(y)
0294 checkedbefore(i) = -y(1);
0295 end
0296 if ~isempty(x) && ~isempty(y);
0297 if x(1) < y(1)
0298 checkedbefore(i) = x(1);
0299 else
0300 checkedbefore(i) = -y(1);
0301 end
0302 end
0303 end
0304
0305
0306
0307
0308
0309
0310
0311
0312
0313
0314 offset = size(LPproblem.A, 2);
0315 for i = 1:length(checkedbefore)
0316 if checkedbefore(i) ==0
0317 continue;
0318 end
0319 pretarget = abs(checkedbefore(i));
0320
0321 if checkedbefore(i) > 0
0322 if any(MILPproblem.A(:,offset+pretarget).*MILPproblem.A(:,offset+i))
0323 display('trouble combining vars'),pause;
0324 end
0325 MILPproblem.A(:,offset+pretarget) = MILPproblem.A(:,offset+pretarget) + MILPproblem.A(:,offset+i);
0326 else
0327 MILPproblem.A(:,offset+pretarget) = MILPproblem.A(:,offset+pretarget) - MILPproblem.A(:,offset+i);
0328 MILPproblem.b = MILPproblem.b - MILPproblem.A(:,offset+i);
0329 end
0330 end
0331
0332 markedforlife = true(size(MILPproblem.A,2), 1);
0333 markedforlife(offset + find(checkedbefore > .5)) = false;
0334
0335 MILPproblem.markedforlife = markedforlife;
0336 MILPproblem.A = MILPproblem.A(:,markedforlife);
0337 MILPproblem.c = MILPproblem.c(markedforlife);
0338 MILPproblem.vartype = MILPproblem.vartype(markedforlife);
0339 MILPproblem.lb = MILPproblem.lb(markedforlife);
0340 MILPproblem.ub = MILPproblem.ub(markedforlife);
0341
0342
0343
0344
0345
0346
0347
0348 end