diff --git a/.github/CONTRIBUTING.md b/.github/CONTRIBUTING.md index e96e5fb1..d5f89431 100644 --- a/.github/CONTRIBUTING.md +++ b/.github/CONTRIBUTING.md @@ -5,9 +5,9 @@ Anybody is welcome to contribute to the development of the RAVEN Toolbox, but pl Each function should start with a commented section describing the function and explaining the parameters. Existing functions can clarify what style should be used. ### Bugfixes, New Features and Functions -* For any development, whether bug fixes, introducing new functions or new/updated features for existing functions: make a separate branch from `devel` and name the branch for instance after the function/feature you are fixing/developing. If you work on a fix, start the branch name with `fix/`, if you work on a feature, start the branch name with `feat/`. Examples: `fix/format_reactions` or `feat/new_algorithms`. +* For any development, whether bug fixes, introducing new functions or new/updated features for existing functions: make a separate branch from `develop` and name the branch for instance after the function/feature you are fixing/developing. If you work on a fix, start the branch name with `fix/`, if you work on a feature, start the branch name with `feat/`. Examples: `fix/format_reactions` or `feat/new_algorithms`. * Make commits to this branch while developing. Aim for backward compatibility, and try to avoid very new MATLAB functions when possible, to accommodate users with older MATLAB versions. -* When you are happy with your new function/feature, make a pull request to the `devel` branch. Also, see [Pull request](#pull-request) below. +* When you are happy with your new function/feature, make a pull request to the `develop` branch. Also, see [Pull request](#pull-request) below. ### Semantic Commits Use semantic commit messages to make it easier to show what you are aiming to do: @@ -28,14 +28,14 @@ A more detailed explanation or comments can be left in the commit description. ### External Software Update * Once the newer version for any external software (BLAST+, DIAMOND, HMMER) is available, identify the newest version which is simultaneously available for macOS, Unix/Linux and Windows -* Create a separate branch from the `devel` branch and name it e.g. `chore/updateBinaries` +* Create a separate branch from the `develop` branch and name it e.g. `chore/updateBinaries` * Commit the changes for each program and operating system separately, e.g. `chore: update HMMER (Win) to 3.2.1` * As soon as all binaries for particular program are updated through commits, update the corresponding license file, if available and place it in e.g. `src/external/blast+` directory. * Update the version number in `RAVENdir/src/external/version.txt` file. -* Do some testing to ensure that the new binaries are working correctly. Upon successful tests, create a Pull Request to the `devel` branch. +* Do some testing to ensure that the new binaries are working correctly. Upon successful tests, create a Pull Request to the `develop` branch. ### Pull Requests -* No changes should be directly committed to the `main` or `devel` branches. Commits are made to side branches, after which pull requests are made for merging with `main` or `devel`. +* No changes should be directly committed to the `main` or `develop` branches. Commits are made to side branches, after which pull requests are made for merging with `main` or `develop`. * The person making the pull request and the one accepting the merge _cannot_ be the same person. * Typically, wait a few days before merging, to allow for other developers to inspect the pull request. * A merge with the mainbranch invokes a new release (see [versioning](#versioning)). @@ -43,11 +43,11 @@ A more detailed explanation or comments can be left in the commit description. ### Versioning RAVEN Toolbox follows [semantic versioning](https://semver.org/). After each successful Pull Request (PR) to `main`, the `version.txt` file should be updated in a separate PR. When this PR is merged, a new [release](https://github.com/CompMetEng/RAVEN/releases) should be specified, with a list detailing the most significant changes since the previous release. -Ensure that the documentation is also updated (in `devel`) before creating a PR to `main`. It can be updated by running `updateDocumentation()`. +Ensure that the documentation is also updated (in `develop`) before creating a PR to `main`. It can be updated by running `updateDocumentation()`. ### Categorize Releases * PRs pushed into `main` should be patch releases by default. -* Developers who push PRs to `devel` may propose target release categories (e.g. minor or patch) that can be commented on and discussed by reviewers and other developers. +* Developers who push PRs to `develop` may propose target release categories (e.g. minor or patch) that can be commented on and discussed by reviewers and other developers. * As a result of the above step, minor releases hopefully can be determined. * To help the process of categorizing releases, some conditions are suggested qualification for minor releases: * with database updates (e.g. KEGG) diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml index c280ae19..7b639ab1 100644 --- a/.github/workflows/documentation.yml +++ b/.github/workflows/documentation.yml @@ -18,7 +18,7 @@ name: Documentation on: push: - branches: [develop3, develop, main] + branches: [develop, main] paths: - '**.m' - '.github/workflows/documentation.yml' diff --git a/docs/v3_review.md b/docs/v3_review.md index 03663bb2..38ae4440 100644 --- a/docs/v3_review.md +++ b/docs/v3_review.md @@ -1,7 +1,6 @@ # RAVEN 3 — development review and roadmap -A review of the RAVEN toolbox aimed at the next major release (**v3**, currently on the -`develop3` branch). The brief: make the codebase **leaner** (remove old/unused/low-value +A review of the RAVEN toolbox for RAVEN 3. The brief: make the codebase **leaner** (remove old/unused/low-value functions) while gaining **a few critical new capabilities** that increase value and usability — with an appetite for ambitious additions. @@ -11,7 +10,7 @@ This document is a proposal and a prioritized to-do list, not a record of work a ## Scope and method -- Reviewed against the **current `develop3` tree**, not `develop`. The four review lenses: +- Reviewed against the RAVEN 3 tree. The four review lenses: **(a)** unnecessary functions to remove, **(b)** missing functionality worth adding, **(c)** existing functionality worth expanding, **(d)** functions worth simplifying. - Method: subsystem-by-subsystem audit of every `.m` file outside `software/`, plus @@ -21,7 +20,7 @@ This document is a proposal and a prioritized to-do list, not a record of work a **COBRA Toolbox (MATLAB)** is wanted; **MetaNetX mapping** is explicitly wanted (an old stale branch exists); new file formats are *not* a priority. -### Already done in `develop3` +### Already done in RAVEN 3 So the review does not repeat completed work. v3 already: @@ -141,14 +140,14 @@ Plus a standalone `mapIDsViaMNXref(type, ids, fromDB, toDB)` ID translator. - **No structural matching** — maps only by existing xref IDs and fuzzy names; ignores the **InChIKey** column, which is the modern primary metabolite key. - **Pre-reorg layout** — lives in `external/MetaNetX/`, carries duplicate `struct_conversion/*` - that has since diverged on develop3. + that has since diverged in RAVEN 3. **Refactor plan** (one focused PR; core is M, the polish items are S each): -1. Relocate to the develop3 layout (e.g. `reconstruction/metanetx/` or `annotation/metanetx/`), - drop the duplicated `struct_conversion/*`, rebase on develop3's `convertMiriams`. +1. Relocate to the RAVEN 3 layout (e.g. `reconstruction/metanetx/` or `annotation/metanetx/`), + drop the duplicated `struct_conversion/*`, rebase on RAVEN 3's `convertMiriams`. 2. **Delete the committed data**; host MNXref 4.5 (gzipped TSVs, and ideally a prebuilt `.mat`) as a `toolbox-assets` release asset and add a `downloadMNXref.m` modeled exactly on the - develop3 KEGG-HMM / `downloadRavenBinaries` on-demand pattern (sentinel-file check, cache + RAVEN 3 KEGG-HMM / `downloadRavenBinaries` on-demand pattern (sentinel-file check, cache under a gitignored dir). 3. Update parsing to the documented **4.5** columns (header-driven `readtable`, handle the `#` comment block and deprecations). diff --git a/src/biomass/scaleBiomassFraction.m b/src/biomass/scaleBiomassFraction.m index 6c66f864..8881f07d 100644 --- a/src/biomass/scaleBiomassFraction.m +++ b/src/biomass/scaleBiomassFraction.m @@ -2,8 +2,8 @@ % scaleBiomassFraction Rescale a biomass component to a target value. % % Rescale a biomass component to a target g/gDW value, optionally -% rescaling a second component so that the biomass total becomes 1 g/gDW, -% whatever the total was before. +% rescaling a second component so that the biomass total becomes +% targetTotal (1 g/gDW by default), whatever the total was before. % The MATLAB counterpart of yeast-GEM's legacy scaleBioMass. % % Parameters @@ -20,10 +20,13 @@ % Name-Value Arguments % -------------------- % balanceOut : char -% Second component name to adjust so the biomass total becomes 1 -% g/gDW. It must be large enough to absorb the change: an error is -% raised if its fraction would become zero or negative, or if it is +% Second component name to adjust so the biomass total becomes +% targetTotal. It must be large enough to absorb the change: an error +% is raised if its fraction would become zero or negative, or if it is % 0 to begin with. Empty / omit to skip balancing. +% targetTotal : double +% Biomass total in g/gDW that balanceOut brings the components to +% (default 1). % % Returns % ------- @@ -33,13 +36,16 @@ % Examples % -------- % model = scaleBiomassFraction(model, biomassConfig, 'protein', 0.5, 'carbohydrate'); +% model = scaleBiomassFraction(model, biomassConfig, 'protein', 0.5, ... +% 'balanceOut', 'carbohydrate', 'targetTotal', 0.95); % % See also % -------- % getBiomassFractions -p=parseRAVENargs(varargin, {'balanceOut',''}); +p=parseRAVENargs(varargin, {'balanceOut',''; 'targetTotal',1}); balanceOut=p.balanceOut; +targetTotal=p.targetTotal; fractions = getBiomassFractions(model, biomassConfig); if ~isfield(fractions, componentName) @@ -68,7 +74,7 @@ ['Cannot balance with %s: current fraction is 0 ' ... '(pseudoreaction missing or empty).'], balanceOut); end - balanceNew = balanceCurrent + (1 - total); + balanceNew = balanceCurrent + (targetTotal - total); if balanceNew <= 0 error('scaleBiomassFraction:negativeBalance', ... ['Cannot balance with %s: its fraction would become %g g/gDW. ' ... diff --git a/src/deprecated/consumeSomething.m b/src/deprecated/consumeSomething.m index 536f8a70..0b468286 100644 --- a/src/deprecated/consumeSomething.m +++ b/src/deprecated/consumeSomething.m @@ -11,9 +11,7 @@ % % findLeakMetabolite inserts allowExcretion (a produce-only option) in % position 4, so the arguments are translated rather than forwarded -% verbatim. Forwarding them verbatim, as this wrapper previously did, -% silently shifted params into allowExcretion and ignoreIntBounds into -% params, so a requested ignoreIntBounds was dropped. +% verbatim. % % See Also % -------- diff --git a/src/gapfilling/gapReport.m b/src/gapfilling/gapReport.m index b28bdc38..7212a86d 100755 --- a/src/gapfilling/gapReport.m +++ b/src/gapfilling/gapReport.m @@ -6,8 +6,8 @@ % Runs a series of connectivity checks on the model - which reactions are % blocked, how many isolated sub-networks exist, which metabolites cannot % be produced (and the minimal set of metabolites that would need to be -% connected to fix that), and, for models with exchange reactions, which -% metabolites can be produced/consumed without any input/output. If +% connected to fix that), and which metabolites can be produced/consumed +% without any input/output, with every boundary reaction closed. If % templateModels is supplied, it also attempts to gap-fill using those % models. A formatted report is printed to the command window; the same % results are returned as outputs for further use. @@ -43,11 +43,11 @@ % matrix where n x m is true if metabolite n allows for production of % metabolite m. % canProduceWithoutInput : cell -% metabolites that could be produced even when there is no input to -% the model. +% metabolites that could be produced with every boundary reaction +% closed. % canConsumeWithoutOutput : cell -% metabolites that could be consumed even when there is no output from -% the model. +% metabolites that could be consumed with every boundary reaction +% closed. % connectedFromTemplates : cell % the reactions that could be connected using the template models. % addedFromTemplates : struct @@ -56,11 +56,15 @@ % % Notes % ----- -% The model must be feasible under its own bounds, and, when the mass -% balancing section runs, also with its boundary closed. Otherwise -% gapReport stops with a RAVEN:infeasibleModel error that names the -% reactions with a forced flux; set the bounds of ATP maintenance and any -% other forced flux to zero first. +% The boundary reactions are those getExchangeRxns returns: exchanges, +% sinks and demands. For a model that keeps its boundary metabolites (an +% unconstrained field), those metabolites are deleted for the connectivity +% checks and kept, balanced, for the closed model. +% +% The model must be feasible under its own bounds and with its boundary +% closed. Otherwise gapReport stops with a RAVEN:infeasibleModel error that +% names the reactions with a forced flux; set the bounds of ATP maintenance +% and any other forced flux to zero first. % % Examples % -------- @@ -84,19 +88,17 @@ fprintf(['Gap analysis for ' esc(model.id) ' - ' esc(model.name) '\n\n']); if isfield(model,'unconstrained') - calculateINOUT=true; closedModel=model; model=simplifyModel(model); else - canConsumeWithoutOutput={}; - canProduceWithoutInput={}; - calculateINOUT=false; + closedModel=model; + [~,boundary]=getExchangeRxns(closedModel); + closedModel.lb(boundary)=0; + closedModel.ub(boundary)=0; end requireFeasible(model,'under its own bounds'); -if calculateINOUT - requireFeasible(closedModel,'once its boundary reactions are closed'); -end +requireFeasible(closedModel,'once its boundary reactions are closed'); model2=model; model2.b=[model2.b inf(numel(model2.mets),1)]; @@ -127,14 +129,12 @@ fprintf(['\t' num2str(i) '. ' esc(minToConnect{i}) '\n']); end -if calculateINOUT==true - fprintf('\n***Mass balancing\n'); - produced=canExchange(closedModel,'produce'); - canProduceWithoutInput=closedModel.mets(produced); - consumed=canExchange(closedModel,'consume'); - canConsumeWithoutOutput=closedModel.mets(consumed); - fprintf([num2str(numel(canConsumeWithoutOutput)) ' metabolites could be consumed without any outputs\n' num2str(numel(canProduceWithoutInput)) ' metabolites could be produced without any inputs\n']); -end +fprintf('\n***Mass balancing\n'); +produced=canExchange(closedModel,'produce'); +canProduceWithoutInput=closedModel.mets(produced); +consumed=canExchange(closedModel,'consume'); +canConsumeWithoutOutput=closedModel.mets(consumed); +fprintf([num2str(numel(canConsumeWithoutOutput)) ' metabolites could be consumed without any outputs\n' num2str(numel(canProduceWithoutInput)) ' metabolites could be produced without any inputs\n']); if ~isempty(templateModels) fprintf('\n***Automated gap-filling\n'); diff --git a/src/installation/fetchToolboxAsset.m b/src/installation/fetchToolboxAsset.m index 1942a794..394b611e 100644 --- a/src/installation/fetchToolboxAsset.m +++ b/src/installation/fetchToolboxAsset.m @@ -35,7 +35,7 @@ websave(localPath,['https://github.com/CompMetEng/toolbox-assets/releases/download/' releaseTag '/' assetName]); catch ME if strcmp(ME.identifier,'MATLAB:webservices:HTTP404StatusCodeError') - error('Failed to download %s, the server returned a 404 error, try again later. If the problem persists please report it on the RAVEN GitHub Issues page: https://github.com/SysBioChalmers/RAVEN/issues',assetName) + error('Failed to download %s, the server returned a 404 error, try again later. If the problem persists please report it on the RAVEN GitHub Issues page: https://github.com/CompMetEng/RAVEN/issues',assetName) else rethrow(ME); end diff --git a/src/solver/getMinNrFluxes.m b/src/solver/getMinNrFluxes.m index 700b4742..655e7287 100755 --- a/src/solver/getMinNrFluxes.m +++ b/src/solver/getMinNrFluxes.m @@ -264,12 +264,6 @@ %One binary per reaction, on the reversible model. Trades a larger %constraint matrix for half the binaries on reversible reactions. -% glpk solver as implemented by COBRA does not work well for MILP. -global CBT_MILP_SOLVER -if strcmp(getpref('RAVEN','solver'),'cobra') && strcmp(CBT_MILP_SOLVER,'glpk') - error('RAVEN:badInput', '%s', 'The current solver is set to ''cobra'', while in COBRA the MILP solver has been set to ''glpk''. The COBRA implementation of glpk is not well suitable for solving MILPs. Please install the Gurobi or an alternative MILP solver.'); -end - exitFlag=1; %The trick to make this possible using a reversible model is that for diff --git a/src/solver/optimizeProb.m b/src/solver/optimizeProb.m index 31991dad..6fc3585e 100755 --- a/src/solver/optimizeProb.m +++ b/src/solver/optimizeProb.m @@ -29,12 +29,11 @@ % % Notes % ----- -% With 'glpk' selected, MILPs are solved by RAVEN's bundled GLPK, and a -% RAVEN:glpkMILP warning recommends Gurobi, which is much faster and more -% reliable on larger MILPs (such as ftINIT or getMinNrFluxes on +% MILPs can be solved with GLPK, either RAVEN's bundled GLPK ('glpk') or +% the COBRA Toolbox's GLPK interface ('cobra' with GLPK as its MILP solver). +% A RAVEN:glpkMILP warning then recommends Gurobi, which is much faster and +% more reliable on larger MILPs (such as ftINIT or getMinNrFluxes on % genome-scale models). -% The GLPK interface of the COBRA Toolbox is not used for MILPs: with -% 'cobra' selected, set a different MILP solver with changeCobraSolver. % % See also % -------- @@ -62,17 +61,11 @@ solver=RAVENSOLVER; if ~all(lower(prob.vartype) == 'c') milp=true; - switch solver - case 'glpk' - warning('RAVEN:glpkMILP', ['Solving a MILP with GLPK. This works for small ' ... - 'problems, but for larger problems Gurobi is much faster and more reliable; ' ... - 'select it with setRavenSolver(''gurobi''). Silence this warning with ' ... - 'warning(''off'',''RAVEN:glpkMILP'').']) - case 'cobra' - if strcmp(CBT_MILP_SOLVER,'glpk') - error(['The GLPK interface of the COBRA Toolbox is not suitable for solving MILPs, ' ... - 'select a different solver with changeCobraSolver() or setRavenSolver().']) - end + if strcmp(solver,'glpk') || (strcmp(solver,'cobra') && strcmp(CBT_MILP_SOLVER,'glpk')) + warning('RAVEN:glpkMILP', ['Solving a MILP with GLPK. This works for small ' ... + 'problems, but for larger problems Gurobi is much faster and more reliable; ' ... + 'select it with setRavenSolver(''gurobi''). Silence this warning with ' ... + 'warning(''off'',''RAVEN:glpkMILP'').']) end else milp=false; diff --git a/tests/function_tests/RavenTestCase.m b/tests/function_tests/RavenTestCase.m index 245e896d..6a5f18a9 100644 --- a/tests/function_tests/RavenTestCase.m +++ b/tests/function_tests/RavenTestCase.m @@ -71,14 +71,18 @@ function assumeDependency(testCase, isAvailable, label) end function assumeMILPSolver(testCase) - % Switch to Gurobi or SCIP for this test, or skip if neither is - % available. + % Switch to Gurobi or SCIP for this test, or else stay on the + % bundled GLPK with its MILP warning silenced. if ~isempty(which('gurobi')) evalc('setRavenSolver(''gurobi'')'); else ok = false; try, evalc('setRavenSolver(''scip'')'); ok = true; catch, end %#ok - testCase.assumeTrue(ok, 'No MILP solver (gurobi/scip) available; test skipped.'); + if ~ok + evalc('setRavenSolver(''glpk'')'); + warnState = warning('off','RAVEN:glpkMILP'); + testCase.addTeardown(@() warning(warnState)); + end end testCase.addTeardown(@() evalc('setRavenSolver(''glpk'')')); end diff --git a/tests/function_tests/tBiomass.m b/tests/function_tests/tBiomass.m index 7b7c0faa..74d1154b 100644 --- a/tests/function_tests/tBiomass.m +++ b/tests/function_tests/tBiomass.m @@ -34,6 +34,24 @@ function scaleBiomassFractionBalancesToOne(testCase) testCase.verifyEqual(f.DNA, 0.3, 'AbsTol', 1e-9); end + function scaleBiomassFractionBalancesToTargetTotal(testCase) + [m, cfg] = tBiomass.twoComponentToyModel(); + out = scaleBiomassFraction(m, cfg, 'protein', 0.7, 'balanceOut', 'DNA', ... + 'targetTotal', 0.9); + f = getBiomassFractions(out, cfg); + testCase.verifyEqual(f.protein, 0.7, 'AbsTol', 1e-9); + testCase.verifyEqual(f.DNA, 0.2, 'AbsTol', 1e-9); + testCase.verifyEqual(f.total, 0.9, 'AbsTol', 1e-9); + end + + function scaleBiomassFractionRefusesNegativeBalanceAtTargetTotal(testCase) + % Protein 0.85 balances against DNA for a total of 1, but not + % for a total of 0.8. + [m, cfg] = tBiomass.twoComponentToyModel(); + testCase.verifyError(@() scaleBiomassFraction(m, cfg, 'protein', 0.85, ... + 'balanceOut', 'DNA', 'targetTotal', 0.8), 'scaleBiomassFraction:negativeBalance'); + end + function scaleBiomassFractionRefusesNegativeBalance(testCase) % Balancing protein 1.2 against DNA would need DNA at -0.2, % which would turn the DNA pseudoreaction around. diff --git a/tests/function_tests/tGapfilling.m b/tests/function_tests/tGapfilling.m index d08bb14c..5f191fff 100644 --- a/tests/function_tests/tGapfilling.m +++ b/tests/function_tests/tGapfilling.m @@ -76,7 +76,9 @@ function fillGapsProducesModel(testCase) end function gapReportRuns(testCase) - evalc('noFluxRxns = gapReport(testCase.model);'); + m = testCase.model; + m.lb(m.lb > 0) = 0; % ATP maintenance, which closing the boundary makes infeasible + evalc('noFluxRxns = gapReport(m);'); testCase.verifyClass(noFluxRxns, 'cell'); end @@ -100,11 +102,30 @@ function gapReportInfeasibleOnceClosedErrors(testCase) testCase.verifySubstring(err.message, 'Maint'); end + function gapReportInfeasibleOnceClosedWithoutUnconstrainedErrors(testCase) + % ExA: =>a Maint: a=>b (forced) ExB: b=>. Feasible with the + % exchanges open; Maint cannot run once they are closed. + m = testCase.tinyModel({'ExA';'Maint';'ExB'}, sparse([1 -1 0; 0 1 -1]), [0;1;0], [10;1000;1000]); + m.mets = {'a';'b'}; m.metNames = m.mets; m.metComps = [1;1]; m.b = zeros(2,1); + err = testCase.gapReportError(m); + testCase.verifySubstring(err.message, 'once its boundary reactions are closed'); + testCase.verifySubstring(err.message, 'Maint'); + end + + function gapReportListsLeaksWithoutUnconstrained(testCase) + % R1: a=>b and R2: b=>a+c make c from nothing once ExA is closed. + m = testCase.tinyModel({'ExA';'R1';'R2'}, sparse([1 -1 1; 0 1 -1; 0 0 1]), [0;0;0], [10;1000;1000]); + m.mets = {'a';'b';'c'}; m.metNames = m.mets; m.metComps = [1;1;1]; m.b = zeros(3,1); + evalc('[~,~,~,~,~,~,produced] = gapReport(m);'); + testCase.verifyTrue(ismember('c', produced)); + end + function gapReportKeepsPercentInModelName(testCase) % A model.name containing "%" must survive intact in the % report header, not be truncated by fprintf misreading it as % a format directive. m = testCase.model; + m.lb(m.lb > 0) = 0; m.name = 'ecoli 50% subset'; out = evalc('gapReport(m);'); testCase.verifySubstring(out, 'ecoli 50% subset'); diff --git a/tests/function_tests/tManipulation.m b/tests/function_tests/tManipulation.m index 6b07eda4..67ec0bcd 100644 --- a/tests/function_tests/tManipulation.m +++ b/tests/function_tests/tManipulation.m @@ -179,8 +179,7 @@ function addTransportOnModelWithoutAnnotationField(testCase) end function addTransportAcceptsRowOrientedMetNames(testCase) - % Regression test for SysBioChalmers/RAVEN#722: a metNames cell - % array keeps whatever orientation the caller passed it in, and + % A metNames cell array keeps whatever orientation the caller passed it in, and % a plain {'a','b'} literal is row-oriented. A single-metabolite % cell (as in addTransportAddsRxn above) cannot catch this --- % a 1x1 cell is the same shape whichever way you call it "row" diff --git a/tests/function_tests/tSolver.m b/tests/function_tests/tSolver.m index dbe99a0b..04645f97 100644 --- a/tests/function_tests/tSolver.m +++ b/tests/function_tests/tSolver.m @@ -59,6 +59,33 @@ function optimizeProbSolvesMILPWithGlpk(testCase) testCase.verifyEqual(res.mipgap, 0); end + function optimizeProbSolvesMILPWithCobraGlpk(testCase) + % The same knapsack through the COBRA Toolbox's GLPK interface. + testCase.assumeDependency(~isempty(which('changeCobraSolver')), 'COBRA Toolbox'); + global CBT_MILP_SOLVER %#ok + origMILP = CBT_MILP_SOLVER; + evalc('changeCobraSolver(''glpk'',''MILP'',0);'); + evalc('setRavenSolver(''cobra'')'); + testCase.addTeardown(@() evalc('setRavenSolver(''glpk'')')); + if ~isempty(origMILP) + testCase.addTeardown(@() changeCobraSolver(origMILP,'MILP',0)); + end + prob = struct(); + prob.c = [-5; -4; -3]; + prob.A = sparse([2 3 1]); + prob.a = prob.A; + prob.b = 4; + prob.csense = 'L'; + prob.lb = zeros(3,1); + prob.ub = ones(3,1); + prob.osense = 1; + prob.vartype = 'BBB'; + res = testCase.verifyWarning(@() optimizeProb(prob,[],false), 'RAVEN:glpkMILP'); + testCase.verifyEqual(res.stat, 1); + testCase.verifyEqual(res.obj, -8, 'AbsTol', 1e-6); + testCase.verifyEqual(res.full, [1; 0; 1], 'AbsTol', 1e-6); + end + function getMinimalMediumRunsOnGlpk(testCase) evalc('setRavenSolver(''glpk'')'); warnState = warning('off','RAVEN:glpkMILP');