0001 function [solution,LPProblem]=solveCobraLPCPLEX(LPProblem,printLevel,basisReuse,conflictResolve,contFunctName,minNorm)
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
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
0065
0066
0067
0068
0069
0070
0071
0072
0073
0074
0075
0076
0077
0078
0079
0080
0081
0082
0083
0084
0085
0086
0087
0088
0089
0090
0091
0092
0093 if ~exist('printLevel','var')
0094 printLevel=2;
0095 end
0096 if ~exist('basisReuse','var')
0097 basisReuse=0;
0098 end
0099 if ~exist('conflictResolve','var')
0100 conflictResolve=0;
0101 end
0102 if ~exist('contFunctName','var')
0103 cpxControl=[];
0104 else
0105 if isstruct(contFunctName)
0106 cpxControl=contFunctName;
0107 else
0108 if ~isempty(contFunctName)
0109
0110
0111
0112 cpxControl=eval(contFunctName);
0113 else
0114 cpxControl=[];
0115 end
0116 end
0117 end
0118 if ~exist('minNorm','var')
0119 minNorm=0;
0120 end
0121
0122 if basisReuse
0123 if isfield(LPProblem,'LPBasis')
0124 basis=LPProblem.LPBasis;
0125
0126 cpxControl.ADVIND=1;
0127 else
0128 basis=[];
0129 end
0130 else
0131 basis=[];
0132
0133 cpxControl.ADVIND=0;
0134 end
0135
0136 if ~isfield(LPProblem,'A')
0137 if ~isfield(LPProblem,'S')
0138 error('Equality constraint matrix must either be a field denoted A or S.')
0139 end
0140 LPProblem.A=LPProblem.S;
0141 end
0142
0143 if ~isfield(LPProblem,'csense')
0144 nMet=size(LPProblem.A);
0145 if printLevel>0
0146 fprintf('%s\n','Assuming equality constraints, i.e. S*v=b');
0147 end
0148
0149 LPProblem.csense(1:nMet,1)='E';
0150 end
0151
0152 if ~isfield(LPProblem,'osense')
0153
0154 LPProblem.osense=-1;
0155 if printLevel>0
0156 fprintf('%s\n','Assuming maximisation of objective');
0157 end
0158 end
0159
0160
0161 [c,x_L,x_U,b,csense,osense] = deal(LPProblem.c,LPProblem.lb,LPProblem.ub,LPProblem.b,LPProblem.csense,LPProblem.osense);
0162
0163 c=full(c*osense);
0164
0165
0166 b=full(b);
0167
0168
0169
0170 if conflictResolve
0171 [m_lin,n]=size(LPProblem.A);
0172 m_quad=0;
0173 m_sos=0;
0174 m_log=0;
0175
0176 mode='full';
0177 fprintf('%s\n%s\n','Building Structure for Conflict Resolution...','...this slows CPLEX down so should not be used for repeated LP');
0178 confgrps = cpxBuildConflict(n,m_lin,m_quad,m_sos,m_log,mode);
0179 prefix=pwd;
0180 suffix='LP_CPLEX_conflict_file.txt';
0181 conflictFile=[prefix '\' suffix];
0182 else
0183 confgrps=[]; conflictFile=[];
0184 end
0185
0186
0187
0188 logfile=[];
0189
0190
0191
0192 savefile=[]; savemode=[];
0193
0194
0195
0196
0197
0198
0199 callback=[];
0200
0201
0202 Prob=[];
0203
0204
0205 IntVars=[]; PI=[]; SC=[]; SI=[]; sos1=[]; sos2=[];
0206
0207
0208 if sum(minNorm)~=0
0209 if length(minNorm)==1
0210
0211 F=speye(length(c))*minNorm;
0212 else
0213 if length(minNorm)~=length(c)
0214 error('Either minNorm is a scalar, or is an n x 1 vector')
0215 else
0216
0217 F=spdiags(minNorm,0,length(c),length(c));
0218 end
0219 end
0220 else
0221 F=[];
0222 end
0223
0224 qc=[];
0225
0226
0227
0228
0229 saRequest =[];
0230
0231
0232 xIP=[];
0233
0234
0235
0236 logcon=[];
0237
0238
0239 tic;
0240
0241
0242 ILOGcomplex=1;
0243 tomlab_cplex=1;
0244 if ~isempty(which('cplexlp')) && tomlab_cplex==0
0245 if ILOGcomplex
0246
0247 if ~isempty(csense)
0248
0249 b_L(csense == 'E',1) = b(csense == 'E');
0250 b_U(csense == 'E',1) = b(csense == 'E');
0251 b_L(csense == 'G',1) = b(csense == 'G');
0252 b_U(csense == 'G',1) = Inf;
0253 b_L(csense == 'L',1) = -Inf;
0254 b_U(csense == 'L',1) = b(csense == 'L');
0255 else
0256 b_L = b;
0257 b_U = b;
0258 end
0259
0260
0261
0262 try
0263 ILOGcplex = Cplex('fba');
0264 catch ME
0265 error('CPLEX not installed or licence server not up')
0266 end
0267
0268 ILOGcplex.Model.sense = 'minimize';
0269
0270
0271 ILOGcplex.Model.obj = c;
0272 ILOGcplex.Model.lb = x_L;
0273 ILOGcplex.Model.ub = x_U;
0274 ILOGcplex.Model.A = LPProblem.A;
0275 ILOGcplex.Model.lhs = b_L;
0276 ILOGcplex.Model.rhs = b_U;
0277
0278 if ~isempty(F)
0279
0280 ILOGcplex.Model.Q=F;
0281 end
0282
0283 if ~isempty(cpxControl)
0284 if isfield(cpxControl,'LPMETHOD')
0285
0286 ILOGcplex.Param.lpmethod.Cur=cpxControl.LPMETHOD;
0287 end
0288 end
0289
0290 if printLevel==0
0291 ILOGcplex.DisplayFunc=[];
0292 else
0293
0294 ILOGcplex.Param.barrier.display.Cur = printLevel;
0295 ILOGcplex.Param.simplex.display.Cur = printLevel;
0296 ILOGcplex.Param.sifting.display.Cur = printLevel;
0297 end
0298
0299
0300 ILOGcplex.solve();
0301
0302 solution.obj = osense*ILOGcplex.Solution.objval;
0303 solution.full = ILOGcplex.Solution.x;
0304 solution.rcost = ILOGcplex.Solution.reducedcost;
0305 solution.dual = ILOGcplex.Solution.dual;
0306 solution.nInfeas = NaN;
0307 solution.sumInfeas = NaN;
0308
0309 solution.origStat = ILOGcplex.Solution.status;
0310 solution.solver = ILOGcplex.Solution.method;
0311 solution.time = ILOGcplex.Solution.time;
0312 else
0313 try
0314 ILOGcplex = Cplex('fba');
0315 catch ME
0316 error('CPLEX not installed or licence server not up')
0317 end
0318
0319 options = cplexoptimset;
0320 switch printLevel
0321 case 0
0322 options = cplexoptimset(options,'Display','off');
0323 case 1
0324 options = cplexoptimset(options,'Display','off');
0325 case 1
0326 options = cplexoptimset(options,'Display','off');
0327 case 1
0328 options = cplexoptimset(options,'Display','off');
0329 end
0330
0331 if ~isempty(csense)
0332 if sum(minNorm)~=0
0333 Aineq = [LPProblem.A(csense == 'L',:); - LPProblem.A(csense == 'G',:)];
0334 bineq = [b(csense == 'L',:); - b(csense == 'G',:)];
0335
0336
0337
0338
0339 [x,fval,exitflag,output,lambda] = cplexqp(F,c,Aineq,bineq,LPProblem.A(csense == 'E',:),b(csense == 'E',1),x_L,x_U,[],options);
0340 else
0341 Aineq = [LPProblem.A(csense == 'L',:); - LPProblem.A(csense == 'G',:)];
0342 bineq = [b(csense == 'L',:); - b(csense == 'G',:)];
0343
0344
0345
0346
0347 [x,fval,exitflag,output,lambda] = cplexlp(c,Aineq,bineq,LPProblem.A(csense == 'E',:),b(csense == 'E',1),x_L,x_U,[],options);
0348 end
0349
0350 solution.obj=osense*fval;
0351 solution.full=x;
0352
0353 solution.dual=lambda.eqlin;
0354 else
0355 Aineq=[];
0356 bineq=[];
0357 if sum(minNorm)~=0
0358 [x,fval,exitflag,output,lambda] = cplexqp(F,c,Aineq,bineq,LPProblem.A,b,x_L,x_U,[],options);
0359 else
0360 [x,fval,exitflag,output,lambda] = cplexlp(c,Aineq,bineq,LPProblem.A,b,x_L,x_U,[],options);
0361 end
0362 solution.obj=osense*fval;
0363 solution.full=x;
0364
0365 solution.dual=sparse(size(LPProblem.A,1),1);
0366 solution.dual(csense == 'E')=lambda.eqlin;
0367
0368 solution.dual(csense == 'L')=lambda.ineqlin(1:nnz(csense == 'L'),1);
0369 solution.dual(csense == 'G')=lambda.ineqlin(nnz(csense == 'L')+1:end,1);
0370 end
0371
0372 solution.rcost=lambda.lower-lambda.upper;
0373 solution.nInfeas = [];
0374 solution.sumInfeas = [];
0375 solution.origStat = output.cplexstatus;
0376 end
0377
0378 Inform = solution.origStat;
0379
0380 else
0381
0382 if ~isempty(csense)
0383
0384 b_L(csense == 'E',1) = b(csense == 'E');
0385 b_U(csense == 'E',1) = b(csense == 'E');
0386 b_L(csense == 'G',1) = b(csense == 'G');
0387 b_U(csense == 'G',1) = Inf;
0388 b_L(csense == 'L',1) = -Inf;
0389 b_U(csense == 'L',1) = b(csense == 'L');
0390 else
0391 b_L = b;
0392 b_U = b;
0393 end
0394
0395
0396
0397
0398
0399 [x, slack, v, rc, f_k, ninf, sinf, Inform, basis] = cplex(c, LPProblem.A, x_L, x_U, b_L, b_U, ...
0400 cpxControl, callback, printLevel, Prob, IntVars, PI, SC, SI, ...
0401 sos1, sos2, F, logfile, savefile, savemode, qc, ...
0402 confgrps, conflictFile, saRequest, basis, xIP, logcon);
0403
0404 solution.full=x;
0405
0406 solution.dual=v;
0407
0408 solution.rcost=rc;
0409 if Inform~=1
0410 solution.obj = NaN;
0411 else
0412 if minNorm==0
0413 solution.obj=f_k*osense;
0414 else
0415 solution.obj=c'*x*osense;
0416 end
0417
0418
0419 end
0420 solution.nInfeas = ninf;
0421 solution.sumInfeas = sinf;
0422 solution.origStat = Inform;
0423 end
0424
0425 timeTaken=NaN;
0426
0427 if Inform~=1 && ~isempty(which('cplex'))
0428 if conflictResolve ==1
0429 if isfield(LPProblem,'mets') && isfield(LPProblem,'rxns')
0430
0431
0432 [nMet,nRxn]=size(LPProblem.A);
0433 totAbbr=nMet+nRxn;
0434 conStrFind=cell(nMet+nRxn,1);
0435 conStrReplace=cell(nMet+nRxn,1);
0436
0437 for m=1:nMet
0438 conStrFind{m,1}=['c' int2str(m) ':'];
0439 conStrReplace{m,1}=[LPProblem.mets{m} ': '];
0440 end
0441
0442 for n=1:nRxn
0443 conStrFind{nMet+n,1}=['x' int2str(n) ' '];
0444 conStrReplace{nMet+n,1}=[LPProblem.rxns{n} ' '];
0445 end
0446 fid1 = fopen(suffix);
0447 fid2 = fopen(['COBRA_' suffix], 'w');
0448 while ~feof(fid1)
0449 tline{1}=fgetl(fid1);
0450
0451
0452
0453 for t=1:totAbbr
0454 tline= strrep(tline, conStrFind{t}, conStrReplace{t});
0455 end
0456 fprintf(fid2,'%s\n', tline{1});
0457 end
0458 fclose(fid1);
0459 fclose(fid2);
0460
0461
0462 else
0463 warning('Need reaction and metabolite abbreviations in order to make a readable conflict resolution file');
0464 end
0465 fprintf('%s\n',['Conflict resolution file written to: ' prefix '\COBRA_' suffix]);
0466 fprintf('%s\n%s\n','The Conflict resolution file gives an irreducible infeasible subset ','of constraints which are making this LP Problem infeasible');
0467 else
0468 if printLevel>0
0469 fprintf('%s\n','No conflict resolution file. Perhaps set conflictResolve = 1 next time.');
0470 end
0471 end
0472 solution.solver = 'cplex_direct';
0473 end
0474
0475
0476
0477
0478
0479
0480 if Inform==1
0481 solution.stat = 1;
0482 if printLevel>0
0483
0484 [ExitText,ExitFlag] = cplexStatus(Inform);
0485 solution.ExitText=ExitText;
0486 solution.ExitFlag=ExitFlag;
0487 fprintf('\n%s%g\n',[ExitText ', Objective '], c'*solution.full*osense);
0488 end
0489 else
0490 if Inform==2
0491 solution.stat = 2;
0492
0493 [ExitText,ExitFlag] = cplexStatus(Inform);
0494 solution.ExitText=ExitText;
0495 solution.ExitFlag=ExitFlag;
0496 fprintf('\n%s%g\n',[ExitText ', Objective '], c'*solution.full*osense);
0497 else
0498 if Inform==3
0499 solution.stat = 0;
0500 else
0501
0502 solution.stat = -1;
0503
0504 [ExitText,ExitFlag] = cplexStatus(Inform);
0505 solution.ExitText=ExitText;
0506 solution.ExitFlag=ExitFlag;
0507 fprintf('\n%s%g\n',[ExitText ', Objective '], c'*solution.full*osense);
0508 end
0509 end
0510 end
0511 solution.time = timeTaken;
0512
0513
0514 if basisReuse
0515 LPProblem.LPBasis=basis;
0516 end
0517
0518 if sum(minNorm)~=0
0519 fprintf('%s\n','This objective corresponds to a flux with minimum Euclidean norm.');
0520 fprintf('%s%d%s\n','The weighting for minimising the norm was ',minNorm,'.');
0521 fprintf('%s\n','Check that the objective is the same without minimising the norm.');
0522 end