Skip to content

Commit

Permalink
fix: swap soplex with scip as solver
Browse files Browse the repository at this point in the history
  • Loading branch information
@edkerk
edkerk committed Feb 8, 2024
1 parent 9b55cf8 commit 7c1069e
Show file tree
Hide file tree
Showing 13 changed files with 317 additions and 245 deletions.
2 changes: 2 additions & 0 deletions .gitignore
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,6 +40,8 @@ Temporary Items
#external/kegg/keggRxns.mat
external/kegg/keggModel.mat
#software
software/scip/
!.keep
#doc

# IDE stuff #
Expand Down
4 changes: 2 additions & 2 deletions INIT/scoreComplexModel.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -218,10 +218,10 @@
% each genes from its fold change between the tissue/celltype(s) in question
% and all other celltypes, or the threshold if supplied. This is a lower
% quality data than protein abundance, since gene abundance is an indirect
% estimate of protiein level. These scores are therefore only used for genes
% estimate of protein level. These scores are therefore only used for genes
% for which there is no HPA data available. The fold changes are transformed
% as min(5*log(x),10) for x > 1 and max(5*log(x),-5) for x < 1 in order to
% have negative scores for lower expressed genes and to scale the scrores
% have negative scores for lower expressed genes and to scale the scores
% to have somewhat lower weights than the HPA scores
tempArrayLevels = arrayData.levels;
tempArrayLevels(isnan(tempArrayLevels)) = 0;
Expand Down
4 changes: 2 additions & 2 deletions doc/INIT/scoreComplexModel.html
View file
Original file line number Diff line number Diff line change
Expand Up @@ -325,10 +325,10 @@ <h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" sr
0218 <span class="comment">% each genes from its fold change between the tissue/celltype(s) in question</span>
0219 <span class="comment">% and all other celltypes, or the threshold if supplied. This is a lower</span>
0220 <span class="comment">% quality data than protein abundance, since gene abundance is an indirect</span>
0221 <span class="comment">% estimate of protiein level. These scores are therefore only used for genes</span>
0221 <span class="comment">% estimate of protein level. These scores are therefore only used for genes</span>
0222 <span class="comment">% for which there is no HPA data available. The fold changes are transformed</span>
0223 <span class="comment">% as min(5*log(x),10) for x &gt; 1 and max(5*log(x),-5) for x &lt; 1 in order to</span>
0224 <span class="comment">% have negative scores for lower expressed genes and to scale the scrores</span>
0224 <span class="comment">% have negative scores for lower expressed genes and to scale the scores</span>
0225 <span class="comment">% to have somewhat lower weights than the HPA scores</span>
0226 tempArrayLevels = arrayData.levels;
0227 tempArrayLevels(isnan(tempArrayLevels)) = 0;
Expand Down
4 changes: 2 additions & 2 deletions doc/installation/checkInstallation.html
View file
Original file line number Diff line number Diff line change
Expand Up @@ -269,7 +269,7 @@ <h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" sr
0207 <a href="#_sub4" class="code" title="subfunction printOrange(stringToPrint)">printOrange</a>(<span class="string">'Fail\n'</span>)
0208 <span class="keyword">end</span>
0209
0210 fprintf([<a href="#_sub2" class="code" title="subfunction str = myStr(InputStr,len)">myStr</a>(<span class="string">' &gt; soplex'</span>,40) <span class="string">'%f'</span>])
0210 fprintf([<a href="#_sub2" class="code" title="subfunction str = myStr(InputStr,len)">myStr</a>(<span class="string">' &gt; scip'</span>,40) <span class="string">'%f'</span>])
0211 <span class="keyword">if</span> res(3).Passed == 1
0212 fprintf(<span class="string">'Pass\n'</span>)
0213 <span class="keyword">else</span>
Expand All @@ -285,7 +285,7 @@ <h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" sr
0223 fprintf([<a href="#_sub2" class="code" title="subfunction str = myStr(InputStr,len)">myStr</a>(<span class="string">' &gt; Set RAVEN solver'</span>,40) <span class="string">'%f'</span>])
0224 <span class="keyword">try</span>
0225 oldSolver=getpref(<span class="string">'RAVEN'</span>,<span class="string">'solver'</span>);
0226 solverIdx=find(strcmp(oldSolver,{<span class="string">'glpk'</span>,<span class="string">'gurobi'</span>,<span class="string">'soplex'</span>,<span class="string">'cobra'</span>}));
0226 solverIdx=find(strcmp(oldSolver,{<span class="string">'glpk'</span>,<span class="string">'gurobi'</span>,<span class="string">'scip'</span>,<span class="string">'cobra'</span>}));
0227 <span class="keyword">catch</span>
0228 solverIdx=0;
0229 <span class="keyword">end</span>
Expand Down
179 changes: 94 additions & 85 deletions doc/solver/optimizeProb.html
View file
Original file line number Diff line number Diff line change
Expand Up @@ -226,93 +226,102 @@ <h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" sr
0172 res.dual = -extra.lambda*prob.osense;
0173 res.rcost = -extra.redcosts*prob.osense;
0174 <span class="comment">%% Use SoPlex</span>
0175 <span class="keyword">case</span> <span class="string">'soplex'</span>
0176 solverparams.tmlim = 1;
0177 solverparams.save = 1;
0178 solverparams = <a href="#_sub1" class="code" title="subfunction s_merged=structUpdate(s_old,s_new)">structUpdate</a>(solverparams,params);
0179 prob.csense = <a href="#_sub2" class="code" title="subfunction paramlist = renameparams(paramlist,old,new)">renameparams</a>(prob.csense, {<span class="string">'L'</span>,<span class="string">'G'</span>,<span class="string">'E'</span>}, {<span class="string">'U'</span>,<span class="string">'L'</span>,<span class="string">'S'</span>});
0180
0181 [ravenDir,currDir]=findRAVENroot();
0182 cd(fullfile(ravenDir,<span class="string">'software'</span>,<span class="string">'GLPKmex'</span>))
0183 glpk(prob.c, prob.A, prob.b, prob.lb, prob.ub, prob.csense, prob.vartype, prob.osense, solverparams);
0184 solverparams.tmlim = defaultparams.timeLimit;
0185 cd(fullfile(ravenDir,<span class="string">'software'</span>,<span class="string">'soplex'</span>))
0186 movefile(fullfile(ravenDir,<span class="string">'software'</span>,<span class="string">'GLPKmex'</span>,<span class="string">'outpb.lp'</span>),<span class="string">'outpb.lp'</span>)
0187 [runCheck, cmdOut] = system([<span class="string">'soplex --solvemode=2 -t'</span> num2str(solverparams.tmlim) <span class="string">' -x outpb.lp &gt; result.out'</span>]);
0188
0189 <span class="keyword">if</span> runCheck ~= 0
0190 error(<span class="string">'SoPlex did not run'</span>)
0191 <span class="keyword">end</span>
0192 <span class="keyword">if</span> verLessThan(<span class="string">'matlab'</span>,<span class="string">'9.9'</span>) <span class="comment">%readlines introduced 2020b</span>
0193 fid=fopen(<span class="string">'result.out'</span>);
0194 line_raw=cell(1000000,1);
0195 i=1;
0196 <span class="keyword">while</span> ~feof(fid)
0197 line_raw{i}=fgetl(fid);
0198 i=i+1;
0199 <span class="keyword">end</span>
0200 line_raw(i:end)=[];
0201 line_raw=string(line_raw);
0202 <span class="keyword">else</span>
0203 <span class="comment">%line_raw=readlines('result.out','EmptyLineRule','skip','WhitespaceRule','trim');</span>
0204 line_raw=readlines(<span class="string">'result.out'</span>);
0205 <span class="keyword">end</span>
0206 delete <span class="string">'result.out'</span>
0207 delete <span class="string">'outpb.lp'</span>
0208 cd(currDir)
0209
0210 <span class="keyword">if</span> find(contains(line_raw,<span class="string">'problem is solved [optimal]'</span>),1) &gt; 0
0211 res.full = zeros(length(prob.c),1);
0212 res.stat = 1;
0213 obj_line = line_raw{contains(line_raw,<span class="string">'Objective value'</span>)};
0214 res.obj = str2double(obj_line(strfind(obj_line,<span class="string">':'</span>)+2:end));
0215 a = find(contains(line_raw,<span class="string">'Primal solution (name, value):'</span>));
0216 z = find(contains(line_raw,<span class="string">'All other variables are zero'</span>));
0217 flux = split(line_raw(a+1:z-1));
0218 rxns = str2double(replace(flux(:,1),<span class="string">'x_'</span>,<span class="string">''</span>));
0219 flux = str2double(flux(:,2));
0220 res.full(rxns) = flux;
0221 <span class="keyword">elseif</span> ~isempty(cmdOut)
0222 error([<span class="string">'SoPlex error: '</span> extractBefore(cmdOut,newline)])
0223 <span class="keyword">else</span>
0224 statusLine = contains(line_raw,<span class="string">'SoPlex status'</span>);
0225 disp(line_raw(statusLine));
0226 res.stat = 0;
0227 res.obj = [];
0228 res.full = [];
0229 <span class="keyword">end</span>
0230
0231 <span class="keyword">otherwise</span>
0232 error(<span class="string">'RAVEN solver not defined or unknown. Try using setRavenSolver(''solver'').'</span>);
0233 <span class="keyword">end</span>
0234 <span class="keyword">if</span> res.stat&gt;0
0235 res.full=res.full(1:size(prob.a,2));
0236 <span class="keyword">end</span>
0237 <span class="keyword">end</span>
0238
0239 <a name="_sub1" href="#_subfunctions" class="code">function s_merged=structUpdate(s_old,s_new)</a>
0240 <span class="comment">%Remove overlapping fields from first struct;</span>
0241 <span class="comment">%Obtain all unique names of remaining fields;</span>
0242 <span class="comment">%Merge both structs</span>
0243 s_merged = rmfield(s_old, intersect(fieldnames(s_old), fieldnames(s_new)));
0244 names = [fieldnames(s_merged); fieldnames(s_new)];
0245 s_merged = cell2struct([struct2cell(s_merged); struct2cell(s_new)], names, 1);
0175 <span class="keyword">case</span> {<span class="string">'soplex'</span>,<span class="string">'scip'</span>} <span class="comment">% Old 'soplex' option also allowed</span>
0176 [xopt,fval,exitflag] = scip([], prob.c, prob.A,-prob.b, prob.b, prob.lb, prob.ub, prob.vartype);
0177
0178 <span class="comment">% [x,fval,exitflag,stats] = scip(H, f, A, rl, ru, lb, ub, xtype, sos, qc, nl, x0, opts)</span>
0179 <span class="comment">%</span>
0180 <span class="comment">% Input arguments*:</span>
0181 <span class="comment">% H - quadratic objective matrix (sparse, optional [NOT TRIL / TRIU])</span>
0182 <span class="comment">% f - linear objective vector</span>
0183 <span class="comment">% A - linear constraint matrix (sparse)</span>
0184 <span class="comment">% rl - linear constraint lhs</span>
0185 <span class="comment">% ru - linear constraint rhs</span>
0186 <span class="comment">% lb - decision variable lower bounds</span>
0187 <span class="comment">% ub - decision variable upper bounds</span>
0188 <span class="comment">% xtype - string of variable integrality ('c' continuous, 'i' integer, 'b' binary)</span>
0189 <span class="comment">% sos - SOS structure with fields type, index and weight (see below)</span>
0190 <span class="comment">% qc - Quadratic Constraints structure with fields Q, l, qrl and qru (see below)</span>
0191 <span class="comment">% nl - Nonlinear Objective and Constraints structure (see below)</span>
0192 <span class="comment">% x0 - primal solution</span>
0193 <span class="comment">% opts - solver options (see below)</span>
0194 <span class="comment">%</span>
0195 <span class="comment">% Return arguments:</span>
0196 <span class="comment">% x - solution vector</span>
0197 <span class="comment">% fval - objective value at the solution</span>
0198 <span class="comment">% exitflag - exit status (see below)</span>
0199 <span class="comment">% stats - statistics structure</span>
0200 <span class="comment">%</span>
0201 <span class="comment">% Option Fields (all optional, see also optiset for a list):</span>
0202 <span class="comment">% solverOpts - specific SCIP options (list of pairs of parameter names and values)</span>
0203 <span class="comment">% maxiter - maximum LP solver iterations</span>
0204 <span class="comment">% maxnodes - maximum nodes to explore</span>
0205 <span class="comment">% maxtime - maximum execution time [s]</span>
0206 <span class="comment">% tolrfun - primal feasibility tolerance</span>
0207 <span class="comment">% display - solver display level [0-5]</span>
0208 <span class="comment">% probfile - write problem to given file</span>
0209 <span class="comment">% presolvedfile - write presolved problem to file</span>
0210 <span class="comment">%</span>
0211 <span class="comment">% Return Status:</span>
0212 <span class="comment">% 0 - Unknown</span>
0213 <span class="comment">% 1 - User Interrupted</span>
0214 <span class="comment">% 2 - Node Limit Reached</span>
0215 <span class="comment">% 3 - Total Node Limit Reached</span>
0216 <span class="comment">% 4 - Stall Node Limit Reached</span>
0217 <span class="comment">% 5 - Time Limit Reached</span>
0218 <span class="comment">% 6 - Memory Limit Reached</span>
0219 <span class="comment">% 7 - Gap Limit Reached</span>
0220 <span class="comment">% 8 - Solution Limit Reached</span>
0221 <span class="comment">% 9 - Solution Improvement Limit Reached</span>
0222 <span class="comment">% 10 - Restart Limit Reached</span>
0223 <span class="comment">% 11 - Problem Solved to Optimality</span>
0224 <span class="comment">% 12 - Problem is Infeasible</span>
0225 <span class="comment">% 13 - Problem is Unbounded</span>
0226 <span class="comment">% 14 - Problem is Either Infeasible or Unbounded</span>
0227
0228 res.origStat = exitflag;
0229 res.full = xopt;
0230 res.obj = fval;
0231
0232 <span class="keyword">switch</span> exitflag
0233 <span class="keyword">case</span> 11
0234 res.stat = 1;
0235 <span class="keyword">case</span> [5, 6, 7, 8, 9, 10, 13]
0236 res.stat = 2;
0237 <span class="keyword">otherwise</span>
0238 res.stat = 0;
0239 <span class="keyword">end</span>
0240 <span class="keyword">otherwise</span>
0241 error(<span class="string">'RAVEN solver not defined or unknown. Try using setRavenSolver(''solver'').'</span>);
0242 <span class="keyword">end</span>
0243 <span class="keyword">if</span> res.stat&gt;0
0244 res.full=res.full(1:size(prob.a,2));
0245 <span class="keyword">end</span>
0246 <span class="keyword">end</span>
0247
0248 <a name="_sub2" href="#_subfunctions" class="code">function paramlist = renameparams(paramlist,old,new)</a>
0249 <span class="keyword">if</span> ~iscell(paramlist)
0250 wasNoCell = true;
0251 paramlist={paramlist};
0252 <span class="keyword">else</span>
0253 wasNoCell = false;
0254 <span class="keyword">end</span>
0255 <span class="keyword">for</span> i=1:numel(old)
0256 paramlist = regexprep(paramlist,old{i},new{i});
0257 <span class="keyword">end</span>
0258 <span class="keyword">if</span> wasNoCell
0259 paramlist=paramlist{1};
0260 <span class="keyword">end</span>
0261 <span class="keyword">end</span></pre></div>
0248 <a name="_sub1" href="#_subfunctions" class="code">function s_merged=structUpdate(s_old,s_new)</a>
0249 <span class="comment">%Remove overlapping fields from first struct;</span>
0250 <span class="comment">%Obtain all unique names of remaining fields;</span>
0251 <span class="comment">%Merge both structs</span>
0252 s_merged = rmfield(s_old, intersect(fieldnames(s_old), fieldnames(s_new)));
0253 names = [fieldnames(s_merged); fieldnames(s_new)];
0254 s_merged = cell2struct([struct2cell(s_merged); struct2cell(s_new)], names, 1);
0255 <span class="keyword">end</span>
0256
0257 <a name="_sub2" href="#_subfunctions" class="code">function paramlist = renameparams(paramlist,old,new)</a>
0258 <span class="keyword">if</span> ~iscell(paramlist)
0259 wasNoCell = true;
0260 paramlist={paramlist};
0261 <span class="keyword">else</span>
0262 wasNoCell = false;
0263 <span class="keyword">end</span>
0264 <span class="keyword">for</span> i=1:numel(old)
0265 paramlist = regexprep(paramlist,old{i},new{i});
0266 <span class="keyword">end</span>
0267 <span class="keyword">if</span> wasNoCell
0268 paramlist=paramlist{1};
0269 <span class="keyword">end</span>
0270 <span class="keyword">end</span></pre></div>
<hr><address>Generated by <strong><a href="http://www.artefact.tk/software/matlab/m2html/" title="Matlab Documentation in HTML">m2html</a></strong> &copy; 2005</address>
</body>
</html>
Loading

0 comments on commit 7c1069e

Please sign in to comment.