From ae6ac2a3db4a356506fc8a84099272bcac1c8fb4 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Tue, 20 Oct 2026 23:14:15 +0200 Subject: [PATCH 1/8] fix: quote YAML 1.1 booleans such as NO in writeYAMLmodel needsQuote now quotes y, n, yes, no, on and off in any case, which YAML 1.1 loaders such as PyYAML read as booleans. This matches the quoting of write_yaml_model in raven-toolbox (CompMetEng/raven-toolbox#178). #10 --- src/io/writeYAMLmodel.m | 5 ++++- tests/function_tests/tIO.m | 26 ++++++++++++++++++++++++++ 2 files changed, 30 insertions(+), 1 deletion(-) diff --git a/src/io/writeYAMLmodel.m b/src/io/writeYAMLmodel.m index 28879551..98dec6a3 100755 --- a/src/io/writeYAMLmodel.m +++ b/src/io/writeYAMLmodel.m @@ -549,7 +549,9 @@ 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 @@ -557,6 +559,7 @@ function emitScalarCore(fid, ~, value, ~, bareNumeric) 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]+)?$|' ... diff --git a/tests/function_tests/tIO.m b/tests/function_tests/tIO.m index 86a06ae6..a0bbbe3f 100644 --- a/tests/function_tests/tIO.m +++ b/tests/function_tests/tIO.m @@ -207,6 +207,32 @@ function writeReadYAMLRoundTrip(testCase) testCase.verifyEqual(numel(m2.rxns), numel(testCase.model.rxns)); end + function writeYAMLmodelQuotesYAML11Booleans(testCase) + % A YAML 1.1 loader such as PyYAML reads a plain y, n, yes, no, + % on or off in any case as a boolean, so these values are + % double-quoted. Longer words that start with them stay plain. + m = testCase.model; + quoted = {'NO'; 'Y'; 'n'; 'on'; 'yes'; 'off'}; + plain = {'NOS1'; 'onion'; 'off-target'}; + m.metNames(1:6) = quoted; + m.metFormulas(1:6) = quoted; + m.metNames(7:9) = plain; + f = [tempname '.yml']; + testCase.addTeardown(@() delete(f)); + evalc('writeYAMLmodel(m, f);'); + lines = splitlines(fileread(f)); + for i = 1:numel(quoted) + testCase.verifyTrue(any(strcmp(lines, [' - name: "' quoted{i} '"'])), quoted{i}); + testCase.verifyTrue(any(strcmp(lines, [' - formula: "' quoted{i} '"'])), quoted{i}); + end + for i = 1:numel(plain) + testCase.verifyTrue(any(strcmp(lines, [' - name: ' plain{i}])), plain{i}); + end + evalc('m2 = readYAMLmodel(f);'); + testCase.verifyEqual(m2.metNames(1:9), m.metNames(1:9)); + testCase.verifyEqual(m2.metFormulas(1:6), quoted); + end + function readYAMLmodelJoinsContinuationStartingWithDash(testCase) % A writer that wraps long scalars can break a name right % before a dash, so the continuation line starts with "-" From da96091d9215b4de9d84e31050deee9e0ed2706b Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 21 Oct 2026 11:39:42 +0000 Subject: [PATCH 2/8] fix: removeLowScoreGenes keeps the first gene in the rule on a tie --- src/manipulation/removeLowScoreGenes.m | 12 ++++++------ tests/function_tests/tINIT.m | 23 +++++++++++++++++++++++ 2 files changed, 29 insertions(+), 6 deletions(-) diff --git a/src/manipulation/removeLowScoreGenes.m b/src/manipulation/removeLowScoreGenes.m index df2378d8..fb9b8c3d 100644 --- a/src/manipulation/removeLowScoreGenes.m +++ b/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 diff --git a/tests/function_tests/tINIT.m b/tests/function_tests/tINIT.m index 56214098..c7bd276a 100644 --- a/tests/function_tests/tINIT.m +++ b/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(); From 41e9162a2b888ab45f309e39f6ac57f5a5d69c4c Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 21 Oct 2026 11:41:07 +0000 Subject: [PATCH 3/8] fix: parseTaskList reads false/0 flags as false, trims ; lists, ignores header case, keeps empty CHANGED bounds as NaN --- src/gapfilling/fitTasks.m | 11 +++-- src/tasks/checkTasks.m | 7 +++- src/tasks/parseTaskList.m | 79 +++++++++++++++++++++++++++-------- tests/function_tests/tTasks.m | 63 ++++++++++++++++++++++++++++ 4 files changed, 136 insertions(+), 24 deletions(-) diff --git a/src/gapfilling/fitTasks.m b/src/gapfilling/fitTasks.m index 78868d16..d311cbbd 100755 --- a/src/gapfilling/fitTasks.m +++ b/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 diff --git a/src/tasks/checkTasks.m b/src/tasks/checkTasks.m index e9e309bc..2c67013d 100755 --- a/src/tasks/checkTasks.m +++ b/src/tasks/checkTasks.m @@ -244,8 +244,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 diff --git a/src/tasks/parseTaskList.m b/src/tasks/parseTaskList.m index 7a177db9..0c946d51 100755 --- a/src/tasks/parseTaskList.m +++ b/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 diff --git a/tests/function_tests/tTasks.m b/tests/function_tests/tTasks.m index 3b858ecf..0cfc3339 100644 --- a/tests/function_tests/tTasks.m +++ b/tests/function_tests/tTasks.m @@ -26,6 +26,69 @@ function parseTaskListReturnsStruct(testCase) testCase.verifyNotEmpty(ts); end + function parseTaskListReadsFlagsListsAndHeadersLeniently(testCase) + % Flags: "false"/"0"/"no" are false, other text is true. Lists + % separated by ";" are trimmed and empty items dropped. Headers + % are matched ignoring case and whitespace. An empty CHANGED + % bound is NaN (leave the bound unchanged). + rows = { + {'' 'id' 'Description' 'Should Fail' 'IN' ' in lb ' 'IN UB' 'OUT' 'OUT LB' 'OUT UB' 'EQU' 'EQU LB' 'EQU UB' 'changed rxn' 'CHANGED LB' 'CHANGED UB' 'PRINT FLUX' 'COMMENTS'} + {'' 'T1' 'first' 'false' 'A[c]; B[c];' '' '5' '' '' '' '' '' '' 'R1; R2' '' '10' '0' ''} + {'' 'T2' 'second' 'yes' ';C[c]' '2' '' '' '' '' '' '' '' 'R3' '-5' '' 'No' ''} + {'' 'T3' 'third' '0' '' '' '' 'D[c]' '' '' '' '' '' '' '' '' 'TRUE' ''} + }; + nCols = 20; + lines = cell(numel(rows),1); + for k = 1:numel(rows) + r = rows{k}; + r(end+1:nCols) = {''}; + lines{k} = strjoin(r, sprintf('\t')); + end + f = [tempname '.txt']; + fid = fopen(f,'w'); + fprintf(fid,'%s\n',lines{:}); + fclose(fid); + cleanup = onCleanup(@() delete(f)); + + ts = parseTaskList(f); + testCase.verifyEqual(numel(ts), 3); + + testCase.verifyFalse(ts(1).shouldFail); + testCase.verifyFalse(ts(1).printFluxes); + testCase.verifyEqual(ts(1).inputs, {'A[c]';'B[c]'}); + testCase.verifyEqual(ts(1).LBin, [0;0]); + testCase.verifyEqual(ts(1).UBin, [5;5]); + testCase.verifyEqual(ts(1).changed, {'R1';'R2'}); + testCase.verifyTrue(all(isnan(ts(1).LBrxn))); + testCase.verifyEqual(ts(1).UBrxn, [10;10]); + + testCase.verifyTrue(ts(2).shouldFail); + testCase.verifyFalse(ts(2).printFluxes); + testCase.verifyEqual(ts(2).inputs, {'C[c]'}); + testCase.verifyEqual(ts(2).LBin, 2); + testCase.verifyEqual(ts(2).LBrxn, -5); + testCase.verifyTrue(isnan(ts(2).UBrxn)); + + testCase.verifyFalse(ts(3).shouldFail); + testCase.verifyTrue(ts(3).printFluxes); + end + + function checkTasksLeavesNaNChangedBoundsAlone(testCase) + m = testCase.taskTestModel(); + task = testCase.taskTestStruct(); + task.changed = {'R2'}; + task.LBrxn = NaN; + task.UBrxn = 1000; + [~, report] = evalc('checkTasks(m, [], true, false, false, task);'); + testCase.verifyTrue(report.ok == 1); + + % R2 can no longer carry flux from a[s] to a[c] with UB 0, and + % its lower bound (-1000) is kept + task.UBrxn = 0; + [~, report] = evalc('checkTasks(m, [], true, false, false, task);'); + testCase.verifyTrue(report.ok == 0); + end + function parseTaskListMatchesExpected(testCase) % The text and Excel task definitions must parse to identical % structs with exactly the expected fields and values. From c49c7da1d165cda0ef6cdf652f9eb9f1571106c4 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 21 Oct 2026 11:42:03 +0000 Subject: [PATCH 4/8] feat: checkTasks warns when a task metabolite matches several model metabolites --- src/tasks/checkTasks.m | 24 ++++++++++++++++++++++++ tests/function_tests/tTasks.m | 10 ++++++++++ 2 files changed, 34 insertions(+) diff --git a/src/tasks/checkTasks.m b/src/tasks/checkTasks.m index 2c67013d..778f4eff 100755 --- a/src/tasks/checkTasks.m +++ b/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 @@ -307,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 diff --git a/tests/function_tests/tTasks.m b/tests/function_tests/tTasks.m index 0cfc3339..2e8feb45 100644 --- a/tests/function_tests/tTasks.m +++ b/tests/function_tests/tTasks.m @@ -26,6 +26,16 @@ function parseTaskListReturnsStruct(testCase) testCase.verifyNotEmpty(ts); end + function checkTasksWarnsWhenSeveralMetabolitesShareAName(testCase) + m = testCase.taskTestModel(); + m.metNames{7} = 'a'; % es is now also a[s], like as + task = testCase.taskTestStruct(); + lastwarn(''); + [~, report] = evalc('checkTasks(m, [], true, false, false, task);'); + testCase.verifySubstring(lastwarn, 'Several model metabolites match'); + testCase.verifyTrue(report.ok == 1); + end + function parseTaskListReadsFlagsListsAndHeadersLeniently(testCase) % Flags: "false"/"0"/"no" are false, other text is true. Lists % separated by ";" are trimmed and empty items dropped. Headers From d43041a5150f2926e61e04811562ea2165e2687c Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 21 Oct 2026 11:42:09 +0000 Subject: [PATCH 5/8] test: use a duplicate that leaves the task feasible --- tests/function_tests/tTasks.m | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tests/function_tests/tTasks.m b/tests/function_tests/tTasks.m index 2e8feb45..11a71eb9 100644 --- a/tests/function_tests/tTasks.m +++ b/tests/function_tests/tTasks.m @@ -28,7 +28,8 @@ function parseTaskListReturnsStruct(testCase) function checkTasksWarnsWhenSeveralMetabolitesShareAName(testCase) m = testCase.taskTestModel(); - m.metNames{7} = 'a'; % es is now also a[s], like as + m.metNames{3} = 'a'; % bc is now also called a[s], like as + m.metComps(3) = 1; task = testCase.taskTestStruct(); lastwarn(''); [~, report] = evalc('checkTasks(m, [], true, false, false, task);'); From feb9a8c052b237469b52c3209a85e3fc52cf723c Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 21 Oct 2026 11:43:07 +0000 Subject: [PATCH 6/8] fix: writeYAMLmodel writes non-integers with the shortest digits that read back exactly --- src/io/writeYAMLmodel.m | 10 +++++++++- tests/function_tests/tIO.m | 22 ++++++++++++++++++++++ 2 files changed, 31 insertions(+), 1 deletion(-) diff --git a/src/io/writeYAMLmodel.m b/src/io/writeYAMLmodel.m index 98dec6a3..b9afd250 100755 --- a/src/io/writeYAMLmodel.m +++ b/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 diff --git a/tests/function_tests/tIO.m b/tests/function_tests/tIO.m index a0bbbe3f..2ca3a5b0 100644 --- a/tests/function_tests/tIO.m +++ b/tests/function_tests/tIO.m @@ -233,6 +233,28 @@ function writeYAMLmodelQuotesYAML11Booleans(testCase) testCase.verifyEqual(m2.metFormulas(1:6), quoted); end + function writeYAMLmodelWritesShortestRoundTripFloats(testCase) + % Non-integers are written with as many digits as needed to read + % back exactly (as Python does), not cut at 15 digits. + m = testCase.model; + m.lb(1) = 1/3; + m.ub(1) = 0.1 + 0.2; + m.lb(2) = 1.5e-7; + m.ub(2) = 2.5; + f = [tempname '.yml']; + testCase.addTeardown(@() delete(f)); + evalc('writeYAMLmodel(m, f);'); + txt = fileread(f); + testCase.verifySubstring(txt, 'lower_bound: 0.3333333333333333'); + testCase.verifySubstring(txt, 'upper_bound: 0.30000000000000004'); + testCase.verifySubstring(txt, 'lower_bound: 1.5e-07'); + testCase.verifySubstring(txt, 'upper_bound: 2.5'); + evalc('m2 = readYAMLmodel(f);'); + testCase.verifyEqual(m2.lb(1), m.lb(1)); + testCase.verifyEqual(m2.ub(1), m.ub(1)); + testCase.verifyEqual(m2.lb(2), m.lb(2)); + end + function readYAMLmodelJoinsContinuationStartingWithDash(testCase) % A writer that wraps long scalars can break a name right % before a dash, so the continuation line starts with "-" From 92a492b28791ab40271dc0cc3af7a6a3da7b2a61 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 21 Oct 2026 11:43:11 +0000 Subject: [PATCH 7/8] test: keep lb below ub in the float round-trip test --- tests/function_tests/tIO.m | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/function_tests/tIO.m b/tests/function_tests/tIO.m index 2ca3a5b0..30f0ff6a 100644 --- a/tests/function_tests/tIO.m +++ b/tests/function_tests/tIO.m @@ -237,16 +237,16 @@ function writeYAMLmodelWritesShortestRoundTripFloats(testCase) % Non-integers are written with as many digits as needed to read % back exactly (as Python does), not cut at 15 digits. m = testCase.model; - m.lb(1) = 1/3; - m.ub(1) = 0.1 + 0.2; + m.lb(1) = 0.1 + 0.2; + m.ub(1) = 1/3; m.lb(2) = 1.5e-7; m.ub(2) = 2.5; f = [tempname '.yml']; testCase.addTeardown(@() delete(f)); evalc('writeYAMLmodel(m, f);'); txt = fileread(f); - testCase.verifySubstring(txt, 'lower_bound: 0.3333333333333333'); - testCase.verifySubstring(txt, 'upper_bound: 0.30000000000000004'); + testCase.verifySubstring(txt, 'upper_bound: 0.3333333333333333'); + testCase.verifySubstring(txt, 'lower_bound: 0.30000000000000004'); testCase.verifySubstring(txt, 'lower_bound: 1.5e-07'); testCase.verifySubstring(txt, 'upper_bound: 2.5'); evalc('m2 = readYAMLmodel(f);'); From 875e45b48b9631077872d135152b5a452ef94535 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 21 Oct 2026 06:28:12 +0000 Subject: [PATCH 8/8] fix: restrict the contractModel call in replaceMets to the metabolite that is left By name, replaceMets deletes the replacement's row and keeps the renamed replaced metabolite. The ids it passed to the final contractModel call (repMetIds) were captured from the replacement only, so after the deletion none of them was found, repIdx was empty and contractModel(model,'mets',[]) checked every reaction. Duplicate reactions that do not involve the merged metabolite were merged as well. repMetIds now also holds the id of the replaced metabolite, so the restriction follows whichever of the two rows survives, and an empty list no longer reaches contractModel. Tests: unrelated duplicate reactions survive a merge by name and by id, and the new duplicate is still contracted. --- src/manipulation/replaceMets.m | 13 +++++-- tests/function_tests/tManipulation.m | 52 ++++++++++++++++++++++++++++ 2 files changed, 63 insertions(+), 2 deletions(-) diff --git a/src/manipulation/replaceMets.m b/src/manipulation/replaceMets.m index bfc31a9f..bcba7d0f 100755 --- a/src/manipulation/replaceMets.m +++ b/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 diff --git a/tests/function_tests/tManipulation.m b/tests/function_tests/tManipulation.m index 67ec0bcd..77e57dfe 100644 --- a/tests/function_tests/tManipulation.m +++ b/tests/function_tests/tManipulation.m @@ -608,6 +608,39 @@ function replaceMetsByNameKeepsBShapeAndRuns(testCase) testCase.verifyEqual(numel(m2.mets), 2); % oxygen and o2 merged end + function replaceMetsByNameLeavesUnrelatedDuplicatesAlone(testCase) + % Merging "oxygen" into "o2" makes r1 a duplicate of r2, which is + % contracted. unrelated1 and unrelated2 are duplicates of each + % other that do not involve either metabolite; contractModel is + % restricted to the reactions of the merged metabolite, so they + % must stay. (By name the replacement's row is the one deleted, + % so the restriction has to follow the metabolite that is left.) + m = tManipulation.duplicateModel(); + evalc('m2 = replaceMets(m, ''oxygen'', ''o2'');'); + testCase.verifyEqual(numel(m2.mets), 3); + testCase.verifyTrue(all(ismember({'unrelated1';'unrelated2'}, m2.rxns))); + testCase.verifyEqual(sum(ismember({'r1';'r2'}, m2.rxns)), 1); + testCase.verifyEqual(numel(m2.rxns), 3); + end + + function replaceMetsByIdLeavesUnrelatedDuplicatesAlone(testCase) + m = tManipulation.duplicateModel(); + evalc('m2 = replaceMets(m, ''a'', ''b'', ''identifiers'', true);'); + testCase.verifyEqual(numel(m2.mets), 3); + testCase.verifyTrue(all(ismember({'unrelated1';'unrelated2'}, m2.rxns))); + testCase.verifyEqual(sum(ismember({'r1';'r2'}, m2.rxns)), 1); + end + + function replaceMetsByNameStillContractsTheNewDuplicate(testCase) + % The restriction must not switch the contraction off: r1 and + % r2 become identical once the two metabolites are one. + m = tManipulation.duplicateModel(); + evalc('m2 = replaceMets(m, ''oxygen'', ''o2'');'); + testCase.verifyLessThan(numel(m2.rxns), numel(m.rxns)); + rowO2 = find(strcmp(m2.metNames, 'o2')); + testCase.verifyEqual(numel(rowO2), 1); + end + function setExchangeBoundsFindsAllMultiExchangeMets(testCase) % Every metabolite exchanged by more than one reaction must be % reported, not just whichever one happens to line up between @@ -897,6 +930,25 @@ function expandModelLeavesDnfAlone(testCase) methods (Static, Access = private) + function m = duplicateModel() + % a and b are named oxygen and o2; r1 (a => x) and r2 (b => x) become + % duplicates once they are one metabolite. unrelated1 and unrelated2 + % (both d =>) are duplicates of each other from the start. + m = struct(); + m.id = 'dup'; + m.rxns = {'r1';'r2';'unrelated1';'unrelated2'}; + m.rxnNames = m.rxns; + m.mets = {'a';'b';'x';'d'}; + m.metNames = {'oxygen';'o2';'x';'d'}; + m.metComps = [1;1;1;1]; + m.comps = {'c'}; + m.compNames = {'cytosol'}; + m.S = sparse([-1 0 0 0; 0 -1 0 0; 1 1 0 0; 0 0 -1 -1]); + m.lb = -1000*ones(4,1); m.ub = 1000*ones(4,1); m.rev = ones(4,1); + m.c = zeros(4,1); m.b = zeros(4,1); + m.genes = {}; m.grRules = repmat({''},4,1); m.rxnGeneMat = sparse(4,0); + end + function m = twoMetModel() % Single reaction a -> b in one compartment. m = struct();