0001 function [solution,model]=solveCobraCPLEX(model,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 if ~exist('printLevel','var')
0087 printLevel=0;
0088 end
0089 if ~exist('basisReuse','var')
0090 basisReuse=0;
0091 end
0092 if ~exist('conflictResolve','var')
0093 conflictResolve=0;
0094 end
0095
0096 if ~exist('minNorm','var')
0097 minNorm=0;
0098 end
0099
0100 if basisReuse
0101 if isfield(model,'LPBasis')
0102 basis=model.LPBasis;
0103
0104
0105 else
0106 basis=[];
0107 end
0108 else
0109 basis=[];
0110
0111 cplex.Param.ADVIND=0;
0112 end
0113
0114 if ~isfield(model,'A')
0115 if ~isfield(model,'S')
0116 error('Equality constraint matrix must either be a field denoted A or S.')
0117 end
0118 model.A=model.S;
0119 end
0120
0121 if ~isfield(model,'csense')
0122 nMet=size(model.A);
0123 if printLevel>0
0124 fprintf('%s\n','Assuming equality constraints, i.e. S*v=b');
0125 end
0126
0127 model.csense(1:nMet,1)='E';
0128 end
0129
0130 if ~isfield(model,'osense')
0131
0132 model.osense=-1;
0133 if printLevel>0
0134 fprintf('%s\n','Assuming maximisation of objective');
0135 end
0136 end
0137
0138 if size(model.A,2)~=length(model.c)
0139 error('dimensions of A & c are inconsistent');
0140 end
0141
0142 if size(model.A,2)~=length(model.lb) || size(model.A,2)~=length(model.ub)
0143 error('dimensions of A & bounds are inconsistent');
0144 end
0145
0146
0147
0148
0149 if conflictResolve
0150 [m_lin,n]=size(model.A);
0151 m_quad=0;
0152 m_sos=0;
0153 m_log=0;
0154
0155 mode='full';
0156 fprintf('%s\n%s\n','Building Structure for Conflict Resolution...','...this slows CPLEX down so should not be used for repeated LP');
0157 confgrps = cpxBuildConflict(n,m_lin,m_quad,m_sos,m_log,mode);
0158 prefix=pwd;
0159 suffix='CPLEX_conflict_file.txt';
0160 conflictFile=[prefix '\' suffix];
0161 else
0162 confgrps=[]; conflictFile=[];
0163 end
0164
0165
0166 try
0167 cplex = Cplex('fba');
0168 catch ME
0169 error('CPLEX not installed or licence server not up')
0170 end
0171
0172
0173 cplex.Model.sense = 'minimize';
0174 if model.osense==1
0175
0176 cplex.Model.obj = model.c;
0177 else
0178
0179 cplex.Model.obj = - model.c;
0180 end
0181
0182 cplex.Model.lb = model.lb;
0183 cplex.Model.ub = model.ub;
0184 cplex.Model.A = model.A;
0185
0186
0187 if isfield(model,'csense')
0188
0189 b_L(model.csense == 'E',1) = model.b(model.csense == 'E');
0190 b_U(model.csense == 'E',1) = model.b(model.csense == 'E');
0191 b_L(model.csense == 'G',1) = model.b(model.csense == 'G');
0192 b_U(model.csense == 'G',1) = Inf;
0193 b_L(model.csense == 'L',1) = -Inf;
0194 b_U(model.csense == 'L',1) = model.b(model.csense == 'L');
0195 cplex.Model.lhs = b_L;
0196 cplex.Model.rhs = b_U;
0197 else
0198 cplex.Model.lhs = model.b;
0199 cplex.Model.rhs = model.b;
0200 end
0201
0202
0203 if sum(minNorm)~=0
0204 if length(minNorm)==1
0205
0206 cplex.Model.Q=model.osense*speye(length(model.c))*minNorm;
0207 else
0208 if length(minNorm)~=length(model.c)
0209 error('Either minNorm is a scalar, or is an n x 1 vector')
0210 else
0211
0212 cplex.Model.Q=model.osense*spdiags(minNorm,0,length(model.c),length(model.c));
0213 end
0214 end
0215 end
0216
0217
0218 if exist('contFunctName','var')
0219 if isstruct(contFunctName)
0220 cplex.Param=contFunctName;
0221 else
0222 if ~isempty(contFunctName)
0223
0224
0225
0226
0227 cplex.Param=Param;
0228 end
0229 end
0230 end
0231
0232 if printLevel==0
0233 cplex.DisplayFunc=[];
0234 else
0235
0236 cplex.Param.barrier.display.Cur = printLevel;
0237 cplex.Param.simplex.display.Cur = printLevel;
0238 cplex.Param.sifting.display.Cur = printLevel;
0239 end
0240
0241
0242 cplex.Param.threads.Cur = 3;
0243
0244
0245 cplex.solve();
0246 solution.origStat = cplex.Solution.status;
0247
0248 if printLevel>0 && solution.origStat~=1
0249
0250 [ExitText,ExitFlag] = cplexStatus(solution.origStat);
0251 solution.ExitText=ExitText;
0252 solution.ExitFlag=ExitFlag;
0253 if any(model.c~=0)
0254 fprintf('\n%s%g\n',[ExitText ', Objective '], model.c'*cplex.Solution.x);
0255 end
0256 end
0257
0258 if solution.origStat==1
0259
0260 solution.obj = model.osense*cplex.Solution.objval;
0261 solution.full = cplex.Solution.x;
0262 solution.rcost = cplex.Solution.reducedcost;
0263 solution.dual = cplex.Solution.dual;
0264 solution.nInfeas = NaN;
0265 solution.sumInfeas = NaN;
0266 solution.solver = cplex.Solution.method;
0267 solution.time = cplex.Solution.time;
0268 else
0269 solution.time=NaN;
0270
0271 if conflictResolve ==1
0272 Cplex.refineConflict
0273 Cplex.writeConflict(suffix)
0274 if isfield(model,'mets') && isfield(model,'rxns')
0275
0276
0277 [nMet,nRxn]=size(model.A);
0278 totAbbr=nMet+nRxn;
0279 conStrFind=cell(nMet+nRxn,1);
0280 conStrReplace=cell(nMet+nRxn,1);
0281
0282 for m=1:nMet
0283 conStrFind{m,1}=['c' int2str(m) ':'];
0284 conStrReplace{m,1}=[model.mets{m} ': '];
0285 end
0286
0287 for n=1:nRxn
0288 conStrFind{nMet+n,1}=['x' int2str(n) ' '];
0289 conStrReplace{nMet+n,1}=[model.rxns{n} ' '];
0290 end
0291 fid1 = fopen(suffix);
0292 fid2 = fopen(['COBRA_' suffix], 'w');
0293 while ~feof(fid1)
0294 tline{1}=fgetl(fid1);
0295
0296
0297
0298 for t=1:totAbbr
0299 tline= strrep(tline, conStrFind{t}, conStrReplace{t});
0300 end
0301 fprintf(fid2,'%s\n', tline{1});
0302 end
0303 fclose(fid1);
0304 fclose(fid2);
0305
0306
0307 else
0308 warning('Need reaction and metabolite abbreviations in order to make a readable conflict resolution file');
0309 end
0310 fprintf('%s\n',['Conflict resolution file written to: ' prefix '\COBRA_' suffix]);
0311 fprintf('%s\n%s\n','The Conflict resolution file gives an irreducible infeasible subset ','of constraints which are making this LP Problem infeasible');
0312 else
0313 if printLevel>0
0314 fprintf('%s\n','No conflict resolution file. Perhaps set conflictResolve = 1 next time.');
0315 end
0316 end
0317 end
0318
0319
0320
0321
0322
0323
0324 if solution.origStat==1
0325 solution.stat = 1;
0326 else
0327
0328 [ExitText,ExitFlag] = cplexStatus(solution.origStat);
0329 solution.ExitText=ExitText;
0330 solution.ExitFlag=ExitFlag;
0331 if any(model.c~=0) && isfield(cplex.Solution,'x')
0332 fprintf('\n%s%g\n',[ExitText ', Objective '], model.c'*cplex.Solution.x);
0333 end
0334 if solution.origStat==2
0335 solution.stat = 2;
0336 else
0337 if solution.origStat==3
0338 solution.stat = 0;
0339 else
0340
0341 solution.stat = -1;
0342 end
0343 end
0344 end
0345
0346
0347 if basisReuse
0348 model.LPBasis=basis;
0349 end
0350
0351 if sum(minNorm)~=0 && printLevel>0
0352 fprintf('%s\n','This objective corresponds to a flux with minimum Euclidean norm.');
0353 if length(minNorm)==1
0354 fprintf('%s%d%s\n','The weighting for minimising the norm was ',minNorm,'.');
0355 else
0356 fprintf('%s%d%s\n','The sum of the weighting for minimising the norm was ',sum(minNorm),'.');
0357 end
0358 fprintf('%s\n','Check that the objective is the same without minimising the norm.');
0359 end