fix: quote YAML 1.1 booleans, fix replaceMets, and align task lists, gene pruning and YAML floats with raven-toolbox

src/gapfilling/fitTasks.m

@@ -262,11 +262,14 @@
         end
         %Add changed bounds
         if ~isempty(taskStructure(i).changed)
-            tModel=setParam(tModel,'lb',taskStructure(i).changed,taskStructure(i).LBrxn);
-            tModel=setParam(tModel,'ub',taskStructure(i).changed,taskStructure(i).UBrxn);
+            %NaN bounds mean "leave unchanged"
+            lbSet=~isnan(taskStructure(i).LBrxn);
+            ubSet=~isnan(taskStructure(i).UBrxn);
+            tModel=setParam(tModel,'lb',taskStructure(i).changed(lbSet),taskStructure(i).LBrxn(lbSet));
+            tModel=setParam(tModel,'ub',taskStructure(i).changed(ubSet),taskStructure(i).UBrxn(ubSet));
             if preMerged
-                tRefModel=setParam(tRefModel,'lb',taskStructure(i).changed,taskStructure(i).LBrxn);
-                tRefModel=setParam(tRefModel,'ub',taskStructure(i).changed,taskStructure(i).UBrxn);
+                tRefModel=setParam(tRefModel,'lb',taskStructure(i).changed(lbSet),taskStructure(i).LBrxn(lbSet));
+                tRefModel=setParam(tRefModel,'ub',taskStructure(i).changed(ubSet),taskStructure(i).UBrxn(ubSet));
             end
         end
 

src/io/writeYAMLmodel.m

@@ -382,7 +382,15 @@ function emitMetaField(fid, key, value)
 if x == floor(x) && abs(x) < 1e16
     s = sprintf('%.1f', x);   % e.g. 1000 -> "1000.0"
 else
-    s = sprintf('%.15g', x);
+    % Shortest text that reads back as exactly x (15 to 17 significant
+    % digits), as Python's repr does, so 1/3 is written 0.3333333333333333
+    % and a value is never rounded on the way out.
+    for digits = 15:17
+        s = sprintf(['%.' num2str(digits) 'g'], x);
+        if str2double(s) == x
+            break
+        end
+    end
 end
 end
 
@@ -549,14 +557,17 @@ function emitScalarCore(fid, ~, value, ~, bareNumeric)
 % because the text itself is not valid as a plain block scalar (leading
 % indicator character, a mid-string ": "/" #", leading/trailing
 % whitespace, or embedded control characters). Implements the YAML 1.2
-% rules for resolving implicit scalar types and analyzing scalar style.
+% rules for resolving implicit scalar types and analyzing scalar style,
+% and also quotes the YAML 1.1 booleans y, n, yes, no, on and off in any
+% case, which YAML 1.1 loaders such as PyYAML read as booleans.
 if isempty(s)
     tf = true;
     return
 end
 
 resolverPattern = [ ...
     '^(true|True|TRUE|false|False|FALSE)$|' ...
+    '^([Yy]|[Nn]|[Yy][Ee][Ss]|[Nn][Oo]|[Oo][Nn]|[Oo][Ff][Ff])$|' ...
     '^[-+]?[0-9][0-9_]*\.[0-9_]*([eE][-+]?[0-9]+)?$|' ...
     '^[-+]?[0-9][0-9_]*[eE][-+]?[0-9]+$|' ...
     '^[-+]?\.[0-9_]+([eE][-+][0-9]+)?$|' ...

src/manipulation/removeLowScoreGenes.m

@@ -19,8 +19,8 @@
 %     score will try to be removed.
 %
 %     If all genes associated with a reaction have a negative score, then the
-%     least-negative gene will remain; if there is a tie, one will be selected
-%     at random.
+%     least-negative gene will remain; if there is a tie, the gene listed first
+%     in the rule will be kept.
 %
 %     If a negative-scoring gene is a subunit in a complex, it will not be
 %     removed; however, the entire complex may be removed. See the following
@@ -149,7 +149,7 @@
 
 
 % get IDs and indices of genes involved in rule
-ruleGenes = unique(regexp(rule,'[^&|\(\) ]+','match'));
+ruleGenes = unique(regexp(rule,'[^&|\(\) ]+','match'),'stable');
 [~,geneInd] = ismember(ruleGenes,genes);
 
 % rules with one or no genes remain unchanged
@@ -164,8 +164,8 @@
     scoreMethod = isozymeScoring;
     negInd = gScores(geneInd) < 0;  % NaNs will return false here
     if all(negInd)
-        % deterministic least-negative: sort descending; stable sort breaks
-        % ties by original gene order (left-to-right in the rule string)
+        % deterministic least-negative: sort descending; the stable sort
+        % breaks ties by the order of the genes in the rule string
         [~,sortOrder] = sort(gScores(geneInd),'descend');
         maxInd = sortOrder(1);
         updatedRule = ruleGenes{maxInd};
@@ -179,7 +179,7 @@
     end
 
     % update ruleGenes and their indices
-    ruleGenes = unique(regexp(updatedRule,'[^&|\(\) ]+','match'));
+    ruleGenes = unique(regexp(updatedRule,'[^&|\(\) ]+','match'),'stable');
     [~,geneInd] = ismember(ruleGenes,genes);
 
 elseif ~contains(rule,'|')  % rule contains enzyme complex

src/manipulation/replaceMets.m

@@ -79,6 +79,12 @@
     error('The to-be-replaced metabolite cannot be found in the model.');
 end
 
+% Either of the two rows can be the one that survives the merge: with identifiers the
+% replacement stays, by name the replaced metabolite is renamed and stays while the
+% replacement's row is deleted below. Both ids are kept so that contractModel is
+% restricted to the metabolite that is left.
+repMetIds = [repMetIds; model.mets(metIdx)];
+
 rxnsWithMet = find(model.S(metIdx,:));
 if verbose==true
     fprintf('\n\nThe following reactions contain the to-be-replaced metabolite as reactant:\n')
@@ -173,10 +179,13 @@
 
 % This could now have created duplicate reactions. Contract model.
 % repIdx may point past the end, or at the wrong row, once idxDelete has
-% shifted things above; re-resolve the replacement metabolite by id.
+% shifted things above; re-resolve the metabolite that is left by id. An empty
+% list must not reach contractModel, which would then check every reaction.
 % distReverse is passed by name rather than positionally as [], which
 % parseRAVENargs would take literally instead of falling back to its
 % documented default of true.
 repIdx = find(ismember(model.mets, repMetIds));
-model=contractModel(model,'mets',repIdx);
+if ~isempty(repIdx)
+    model=contractModel(model,'mets',repIdx);
+end
 end

src/tasks/checkTasks.m

@@ -112,6 +112,13 @@
     if ~isempty(taskStructure(i).inputs)
         [I, J]=ismember(upper(taskStructure(i).inputs),modelMets);
         J=J(I); %Only keep the ones with matches
+        %Several model metabolites can share a name[comp]. Only the first
+        %one is constrained, so tell the user when that happens
+        dupMets=duplicateMetNames(taskStructure(i).inputs(I),modelMets);
+        if ~isempty(dupMets)
+            EM=['Several model metabolites match input(s) in "[' taskStructure(i).id '] ' taskStructure(i).description '": ' strjoin(dupMets,', ') '. Only the first match is constrained'];
+            warning('RAVEN:warning', '%s', EM);
+        end
         K=ismember(upper(taskStructure(i).inputs),'ALLMETS');
         L=~cellfun('isempty',strfind(upper(taskStructure(i).inputs),'ALLMETSIN'));
         %Check that all metabolites are either real metabolites or
@@ -165,6 +172,13 @@
     if ~isempty(taskStructure(i).outputs)
         [I, J]=ismember(upper(taskStructure(i).outputs),modelMets);
         J=J(I); %Only keep the ones with matches
+        %Several model metabolites can share a name[comp]. Only the first
+        %one is constrained, so tell the user when that happens
+        dupMets=duplicateMetNames(taskStructure(i).outputs(I),modelMets);
+        if ~isempty(dupMets)
+            EM=['Several model metabolites match output(s) in "[' taskStructure(i).id '] ' taskStructure(i).description '": ' strjoin(dupMets,', ') '. Only the first match is constrained'];
+            warning('RAVEN:warning', '%s', EM);
+        end
         K=ismember(upper(taskStructure(i).outputs),'ALLMETS');
         L=~cellfun('isempty',strfind(upper(taskStructure(i).outputs),'ALLMETSIN'));
         %Check that all metabolites are either real metabolites or
@@ -244,8 +258,11 @@
     end
     %Add changed bounds
     if ~isempty(taskStructure(i).changed)
-        tModel=setParam(tModel,'lb',taskStructure(i).changed,taskStructure(i).LBrxn);
-        tModel=setParam(tModel,'ub',taskStructure(i).changed,taskStructure(i).UBrxn);
+        %NaN bounds mean "leave unchanged"
+        lbSet=~isnan(taskStructure(i).LBrxn);
+        ubSet=~isnan(taskStructure(i).UBrxn);
+        tModel=setParam(tModel,'lb',taskStructure(i).changed(lbSet),taskStructure(i).LBrxn(lbSet));
+        tModel=setParam(tModel,'ub',taskStructure(i).changed(ubSet),taskStructure(i).UBrxn(ubSet));
     end
 
     %Solve and print
@@ -304,3 +321,13 @@
 taskReport.ok=ok;
 
 end
+
+function dupMets = duplicateMetNames(names,modelMets)
+% Return the (upper-case) names that match more than one model metabolite.
+dupMets = unique(upper(names));
+isDup = false(size(dupMets));
+for k = 1:numel(dupMets)
+    isDup(k) = sum(strcmp(modelMets,dupMets{k})) > 1;
+end
+dupMets = dupMets(isDup);
+end

src/tasks/parseTaskList.m

@@ -41,14 +41,20 @@
 %       Several IDs can be delimited by ";". If so, then the same bounds
 %       are used for all reactions. If that is not wanted, then use several
 %       rows for the task
-%     - CHANGED LB : lower bound for the reaction
-%     - CHANGED UB : upper bound for the reaction
+%     - CHANGED LB : lower bound for the reaction. If empty, the lower bound
+%       of the reaction is left as it is
+%     - CHANGED UB : upper bound for the reaction. If empty, the upper bound
+%       of the reaction is left as it is
 %     - SHOULD FAIL : true if the correct behavior of the model is to not
-%       have a feasible solution given the constraints (optional, default
+%       have a feasible solution given the constraints. Any non-empty value
+%       except 0, false and no (any case) counts as true (optional, default
 %       false)
 %     - PRINT FLUX : true if the function should print the corresponding
-%       flux distribution for a task. Can be useful for testing (optional,
-%       default false)
+%       flux distribution for a task. Can be useful for testing. Read in
+%       the same way as SHOULD FAIL (optional, default false)
+%
+%     Headers are matched ignoring case and surrounding whitespace. In the
+%     ";"-delimited cells, items are trimmed and empty items are dropped.
 %
 % Returns
 % -------
@@ -71,8 +77,10 @@
 %       reversible and 0 for irreversible)
 %     - UBequ : array with upper bounds on equations (default, 1000)
 %     - changed : cell array with reactions to change bounds for
-%     - LBrxn : array with lower bounds on changed reactions
-%     - UBrxn : array with upper bounds on changed reactions
+%     - LBrxn : array with lower bounds on changed reactions (NaN if empty,
+%       meaning the bound is left unchanged)
+%     - UBrxn : array with upper bounds on changed reactions (NaN if empty,
+%       meaning the bound is left unchanged)
 %
 % Examples
 % --------
@@ -136,8 +144,12 @@
 
 %Captions
 columns={'ID';'DESCRIPTION';'IN';'IN LB';'IN UB';'OUT';'OUT LB';'OUT UB';'EQU';'EQU LB';'EQU UB';'CHANGED RXN';'CHANGED LB';'CHANGED UB';'SHOULD FAIL';'PRINT FLUX';'COMMENTS'};
-%Match the columns
-[I, colI]=ismember(columns,raw(1,:));
+%Match the columns. Headers are matched ignoring case and surrounding
+%whitespace, so "id" and " IN LB " are accepted.
+header=raw(1,:);
+isText=cellfun(@ischar,header);
+header(isText)=upper(strtrim(header(isText)));
+[I, colI]=ismember(columns,header);
 
 %Only ID is required. Columns that are absent get index 0 from ismember, which
 %cannot be used to index raw, so point them at a spare empty column instead:
@@ -204,10 +216,10 @@
     task.id=raw{2,colI(1)};
 end
 task.description=raw{2,colI(2)};
-if ~isempty(raw{2,colI(15)})
+if isTruthy(raw{2,colI(15)})
     task.shouldFail=true;
 end
-if ~isempty(raw{2,colI(16)})
+if isTruthy(raw{2,colI(16)})
     task.printFluxes=true;
 end
 if ~isempty(raw{2,colI(17)})
@@ -217,14 +229,14 @@
 for i=2:size(raw,1)
     %Set the inputs
     if ischar(raw{i,colI(3)})
-        inputs=regexp(raw{i,colI(3)},';','split');
+        inputs=splitList(raw{i,colI(3)});
         task.inputs=[task.inputs;inputs(:)];
         task.LBin=[task.LBin;ones(numel(inputs),1)*raw{i,colI(4)}];
         task.UBin=[task.UBin;ones(numel(inputs),1)*raw{i,colI(5)}];
     end
     %Set the outputs
     if ischar(raw{i,colI(6)})
-        outputs=regexp(raw{i,colI(6)},';','split');
+        outputs=splitList(raw{i,colI(6)});
         task.outputs=[task.outputs;outputs(:)];
         task.LBout=[task.LBout;ones(numel(outputs),1)*raw{i,colI(7)}];
         task.UBout=[task.UBout;ones(numel(outputs),1)*raw{i,colI(8)}];
@@ -249,10 +261,11 @@
     end
     %Add changed bounds
     if ischar(raw{i,colI(12)})
-        changed=regexp(raw{i,colI(12)},';','split');
+        changed=splitList(raw{i,colI(12)});
         task.changed=[task.changed;changed(:)];
-        task.LBrxn=[task.LBrxn;ones(numel(changed),1)*raw{i,colI(13)}];
-        task.UBrxn=[task.UBrxn;ones(numel(changed),1)*raw{i,colI(14)}];
+        %An empty bound means "leave this bound as it is", stored as NaN
+        task.LBrxn=[task.LBrxn;ones(numel(changed),1)*boundOrNaN(raw{i,colI(13)})];
+        task.UBrxn=[task.UBrxn;ones(numel(changed),1)*boundOrNaN(raw{i,colI(14)})];
     end
 
     %Check if it should add more constraints
@@ -271,10 +284,10 @@
             task.id=raw{i+1,colI(1)};
         end
         task.description=raw{i+1,colI(2)};
-        if ~isempty(raw{i+1,colI(15)})
+        if isTruthy(raw{i+1,colI(15)})
             task.shouldFail=true;
         end
-        if ~isempty(raw{i+1,colI(16)})
+        if isTruthy(raw{i+1,colI(16)})
             task.printFluxes=true;
         end
         if ~isempty(raw{i+1,colI(17)})
@@ -286,3 +299,33 @@
 %Should add more checks, such as unique IDs and missing headers
 
 end
+
+function tf = isTruthy(x)
+% A flag cell is true unless it is empty, 0, or one of the texts "0", "false"
+% and "no" (any case).
+if isempty(x)
+    tf = false;
+elseif ischar(x) || isstring(x)
+    tf = ~ismember(lower(strtrim(char(x))),{'','0','false','no'});
+elseif isnumeric(x) || islogical(x)
+    tf = any(x(:) ~= 0);
+else
+    tf = true;
+end
+end
+
+function items = splitList(str)
+% Split a ";"-delimited cell on ";", trim each item and drop empty ones, so
+% that "A; B;" gives {'A';'B'}.
+items = strtrim(regexp(str,';','split'));
+items = items(~cellfun(@isempty,items));
+end
+
+function b = boundOrNaN(x)
+% Bound for a changed reaction: NaN (leave unchanged) when the cell is empty.
+if isempty(x)
+    b = NaN;
+else
+    b = x;
+end
+end

tests/function_tests/tINIT.m

@@ -62,6 +62,29 @@ function removeLowScoreGenesKeepsFieldsAligned(testCase)
             testCase.verifyEqual(full(mm.rxnGeneMat), ones(1, numel(mm.genes)));
         end
 
+        function removeLowScoreGenesTieKeepsFirstGeneInRule(testCase)
+            % With every isozyme negative and tied, the gene listed first in
+            % the rule is kept (not the alphabetically first one), and the
+            % order of the surviving genes follows the rule.
+            m = struct();
+            m.id = 'test';
+            m.rxns = {'R1';'R2'}; m.rxnNames = {'R1';'R2'};
+            m.mets = {'A';'B'}; m.metNames = {'A';'B'}; m.metComps = [1;1];
+            m.comps = {'c'}; m.compNames = {'c'};
+            m.S = sparse([-1 -1;1 1]); m.lb = [0;0]; m.ub = [1000;1000];
+            m.rev = [0;0]; m.c = [0;0]; m.b = [0;0];
+            m.genes = {'G1';'G2';'G3';'G4'};
+            m.grRules = {'G2 or G1';'G4 or G3 or G2'};
+            m.rxnGeneMat = sparse([1 1 0 0;0 1 1 1]);
+            evalc('[mm, removed] = removeLowScoreGenes(m, [-1;-1;1;-1]);');
+            testCase.verifyEqual(mm.grRules{1}, 'G2');
+            testCase.verifyEqual(mm.grRules{2}, 'G3');
+            testCase.verifyEqual(sort(removed), {'G1';'G4'});
+
+            evalc('mm = removeLowScoreGenes(m, [1;1;1;1]);');
+            testCase.verifyEqual(mm.grRules, m.grRules);
+        end
+
         function reverseRxnsRuns(testCase)
             % R1 = '=> a[s]'; R3 = 'a[c] <=> b[c] + c[c]'
             testModel = getTstModel();