Follow-ups from the RAVEN 3 Guide review

.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,26 +28,26 @@ 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)).
 
 ### 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)

.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'

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).

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. ' ...

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
 % --------

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');

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

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

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;

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<NOCOM>
-                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

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.