From c4e922288111261e89ba85930c176730cf40c263 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 28 Oct 2026 11:17:25 +0100 Subject: [PATCH 01/27] ftINIT: resolveTies takes each step's fluxes from a least-total-flux solution After the id-rank phase, a further phase holds the on/off state and minimises the summed absolute flux, so the direction and magnitude at which the next step forces a reaction do not depend on the solver's vertex. --- CHANGELOG.md | 6 ++++ src/INIT/ftINIT.m | 7 +++-- src/INIT/ftINITInternalAlg.m | 59 ++++++++++++++++++++++++++++++++++-- tests/function_tests/tINIT.m | 31 +++++++++++++++++++ 4 files changed, 99 insertions(+), 4 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 41cdf996..9f83d30d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,11 @@ # Changelog +## Unreleased + +* `ftINIT` with `resolveTies`: each main-extraction step takes its fluxes from the + least-total-flux solution for the tie-broken on/off state, so the direction and magnitude + at which the next step forces a reaction no longer depend on the solver's vertex. + ## 3.0.0b1 — 2026-09-07 First beta toward 3.0. Changes below are relative to RAVEN 2.9.3, the last 2.x release: diff --git a/src/INIT/ftINIT.m b/src/INIT/ftINIT.m index 1c66ba91..b5df2da3 100644 --- a/src/INIT/ftINIT.m +++ b/src/INIT/ftINIT.m @@ -71,8 +71,11 @@ % if true, pin each step's degenerate MILP to a canonical answer instead % of relying on the solver Seed alone: hold the score objective at its % optimum, then minimise the count of removable ("on") reactions, then, -% among the sparsest, minimise their summed reaction-id rank. Applies to -% every main-extraction step (ftINITInternalAlg) and to the final +% among the sparsest, minimise their summed reaction-id rank; in each +% main-extraction step (ftINITInternalAlg) the fluxes, which set the +% direction and magnitude a reaction is forced at in the next step, then +% come from the least-total-flux solution for that on/off state. Applies +% to every main-extraction step and to the final % gap-filling step (fitTasks' resolveTies) -- see raven-gecko-parity#104 % and raven-toolbox's resolve_ties (default false). % proveAbsGap : double diff --git a/src/INIT/ftINITInternalAlg.m b/src/INIT/ftINITInternalAlg.m index 9e8ebc84..13af878a 100644 --- a/src/INIT/ftINITInternalAlg.m +++ b/src/INIT/ftINITInternalAlg.m @@ -44,8 +44,11 @@ % rank. Positive-score reactions carry a continuous, not a % free binary, on/off variable (pinned near 1 by the score % objective), so only the negative-score ones are a genuine -% tie to resolve -- matches raven-toolbox's resolve_ties / -% _resolve_ties (raven-gecko-parity#104) (optional, default +% tie to resolve. The fluxes are then taken from the +% least-total-flux solution for that on/off state, so the +% directions and forced magnitudes the next ftINIT step +% reads from them are pinned too -- matches raven-toolbox's +% resolve_ties (raven-gecko-parity#104) (optional, default % false) % % deletedRxns reactions which were deleted by the algorithm (only @@ -522,6 +525,58 @@ if isFeasible2b %Adopt the tie-resolved solution for the on/off read-out below. res = res2b; + + %Phase 2c: the least total absolute flux for this on/off state. + %The next step forces each turned-on reaction in the direction + %of its flux, at min(0.99*|flux|,0.1), so the flux read-out + %comes from a minimal-flux solution instead of an arbitrary + %vertex. A binary on/off variable is fixed at its value; a + %continuous one that is on keeps at least its value, one that + %is off at most its value. Directions of reversible scored + %reactions stay free, so this is still a MILP. + nCols = numel(prob.c); + minProb = tieProb2b; + if isfield(minProb,'start') + minProb = rmfield(minProb,'start'); + end + x = res2b.full; + negVal = round(x(negCols)); + minProb.lb(negCols) = negVal; + minProb.ub(negCols) = negVal; + posCols = [onoffPosIrrev, onoffPosRev]; + posVal = min(max(x(posCols),0),1); + posOn = posVal >= 0.5; + minProb.lb(posCols(posOn)) = posVal(posOn); + minProb.ub(posCols(~posOn)) = posVal(~posOn); + %|v| is v for a forward-only flux, -v for a backward-only one, + %and an added t >= v, t >= -v otherwise. + fwdOnly = minProb.lb(1:nRxns) >= 0; + bwdOnly = ~fwdOnly & minProb.ub(1:nRxns) <= 0; + tRxns = find(~fwdOnly & ~bwdOnly); + nT = numel(tRxns); + selV = sparse(1:nT, tRxns, 1, nT, nCols); + minProb.A = [minProb.A sparse(size(minProb.A,1),nT); -selV speye(nT); selV speye(nT)]; + minProb.a = minProb.A; + minProb.b = [minProb.b; zeros(2*nT,1)]; + minProb.csense = [minProb.csense repmat('G',1,2*nT)]; + minProb.lb = [minProb.lb; zeros(nT,1)]; + minProb.ub = [minProb.ub; inf(nT,1)]; + minProb.vartype = [minProb.vartype repmat('C',1,nT)]; + minProb.c = [zeros(nCols,1); ones(nT,1)]; + minProb.c(fwdOnly) = 1; + minProb.c(bwdOnly) = -1; + minParams = params; + minParams.MIPGap = 0; + minParams.MIPGapAbs = 1e-6; + res2c = optimizeProb(minProb,minParams,verbose); + [isFeasible2c, isOptimal2c] = checkSolution(res2c); + if ~isOptimal2c + unproven{end+1} = 'phase2c-minflux'; + end + if isFeasible2c + %Only the fluxes change; the status stays that of phase 2b. + res.full = res2c.full(1:nCols); + end end end if ~isempty(unproven) diff --git a/tests/function_tests/tINIT.m b/tests/function_tests/tINIT.m index 4564e29d..cdbd1437 100644 --- a/tests/function_tests/tINIT.m +++ b/tests/function_tests/tINIT.m @@ -384,6 +384,37 @@ function ftINITResolveTiesAndProveAbsGapStillRunCorrectly(testCase) testCase.verifyTrue(all(strcmp(resBoth.rxns, expected))); end + function ftINITInternalAlgResolveTiesGivesMinimalFluxes(testCase) + % With resolveTies the fluxes come from the least-total-flux + % solution for the tie-broken on/off state. The positive + % reversible R3 (m <=> q) can be on in either direction: backward + % (import q, feed the essential E) carries 0.4 in total, forward + % (a via the free R1 into R3 and E) 0.8. The free loop L1/L2 + % carries nothing and the negative-score dead end R2 stays off. + testCase.assumeMILPSolver(); + m.rxns = {'R1';'R2';'R3';'E';'L1';'L2';'EX_a';'EX_p';'EX_q'}; + m.mets = {'a';'m';'p';'q';'n';'z'}; + % R1 R2 R3 E L1 L2 EXa EXp EXq + m.S = sparse([-1 -1 0 0 0 0 -1 0 0; % a + 1 0 -1 -1 -1 1 0 0 0; % m + 0 0 0 1 0 0 0 -1 0; % p + 0 0 1 0 0 0 0 0 -1; % q + 0 0 0 0 1 -1 0 0 0; % n + 0 1 0 0 0 0 0 0 0]); % z + m.lb = [0;0;-1000;0;0;0;-1000;-1000;-1000]; + m.ub = 1000*ones(9,1); + m.rev = double(m.lb < 0); + scores = [0;-1;1;0;0;0;0;0;0]; + fluxIn = 0.1*ones(9,1); + [deleted,~,~,turnedOn,fluxes] = ftINITInternalAlg(m,scores,[],{'E'}, ... + [],false,false,struct(),[],fluxIn,false,0.1,true); + testCase.verifyTrue(turnedOn(3)); + testCase.verifyTrue(deleted(2)); + % E is forced at min(0.99*0.1,0.1) and needs the 0.1 of R3. + expected = [0;0;-0.1;0.1;0;0;0;0.1;-0.1]; + testCase.verifyEqual(fluxes, expected, 'AbsTol', 1e-6); + end + function ftINITWithTaskRuns(testCase) % A task requiring e[s] from a[s] forces R2 and R7 to be essential. testCase.assumeMILPSolver(); From d562aee99a3ea18ce81efea671a76cbdf4d6ac7e Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 28 Oct 2026 18:16:42 +0100 Subject: [PATCH 02/27] ftINIT: the resolveTies minimal-flux phase stops at a 1e-6 relative or absolute gap --- CHANGELOG.md | 1 + src/INIT/ftINITInternalAlg.m | 4 +++- 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 9f83d30d..ecd42a05 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,7 @@ * `ftINIT` with `resolveTies`: each main-extraction step takes its fluxes from the least-total-flux solution for the tie-broken on/off state, so the direction and magnitude at which the next step forces a reaction no longer depend on the solver's vertex. + That phase stops at a relative or absolute gap of 1e-6. ## 3.0.0b1 — 2026-09-07 diff --git a/src/INIT/ftINITInternalAlg.m b/src/INIT/ftINITInternalAlg.m index 13af878a..47fe78fe 100644 --- a/src/INIT/ftINITInternalAlg.m +++ b/src/INIT/ftINITInternalAlg.m @@ -534,6 +534,8 @@ %continuous one that is on keeps at least its value, one that %is off at most its value. Directions of reversible scored %reactions stay free, so this is still a MILP. + %It stops at a relative or absolute gap of 1e-6, as raven-toolbox's + %phase 2c. nCols = numel(prob.c); minProb = tieProb2b; if isfield(minProb,'start') @@ -566,7 +568,7 @@ minProb.c(fwdOnly) = 1; minProb.c(bwdOnly) = -1; minParams = params; - minParams.MIPGap = 0; + minParams.MIPGap = 1e-6; minParams.MIPGapAbs = 1e-6; res2c = optimizeProb(minProb,minParams,verbose); [isFeasible2c, isOptimal2c] = checkSolution(res2c); From 22d24878c6562852be6faf5d030dbd55933aa9a0 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 28 Oct 2026 18:23:46 +0100 Subject: [PATCH 03/27] ftINIT: resolveTies fixes equal-total-flux directions by id rank and ranks a merged reaction by its smallest original id --- CHANGELOG.md | 5 +++ src/INIT/ftINIT.m | 22 ++++++++++--- src/INIT/ftINITInternalAlg.m | 60 +++++++++++++++++++++++++++++++----- tests/function_tests/tINIT.m | 30 ++++++++++++++++++ 4 files changed, 106 insertions(+), 11 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index ecd42a05..2e5acd85 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,11 @@ least-total-flux solution for the tie-broken on/off state, so the direction and magnitude at which the next step forces a reaction no longer depend on the solver's vertex. That phase stops at a relative or absolute gap of 1e-6. +* `ftINIT` with `resolveTies`: among the least-total-flux solutions, each main-extraction + step takes its fluxes from the one with the least id-rank-weighted backward flux, so the + direction of a flux that carries the same total either way no longer depends on the + solver. The id-rank phases rank a merged reaction by the smallest original reaction id in + its group, so they no longer depend on which member names it. ## 3.0.0b1 — 2026-09-07 diff --git a/src/INIT/ftINIT.m b/src/INIT/ftINIT.m index b5df2da3..6efa2e47 100644 --- a/src/INIT/ftINIT.m +++ b/src/INIT/ftINIT.m @@ -74,8 +74,10 @@ % among the sparsest, minimise their summed reaction-id rank; in each % main-extraction step (ftINITInternalAlg) the fluxes, which set the % direction and magnitude a reaction is forced at in the next step, then -% come from the least-total-flux solution for that on/off state. Applies -% to every main-extraction step and to the final +% come from the least-total-flux solution for that on/off state and, +% among those, from the one with the least id-rank-weighted backward +% flux (a merged reaction ranked by the smallest original id in its +% group). Applies to every main-extraction step and to the final % gap-filling step (fitTasks' resolveTies) -- see raven-gecko-parity#104 % and raven-toolbox's resolve_ties (default false). % proveAbsGap : double @@ -194,6 +196,18 @@ %essential rxn. fluxes = ones(length(prepData.minModel.rxns), 1).*0.1; +%resolveTies ranks a merged reaction by its group's smallest original id, +%which does not depend on which member names the merged reaction. +rankKeys = prepData.minModel.rxns; +[inRef, refInd] = ismember(rankKeys, prepData.refModel.rxns); +for i = find(inRef)' + g = prepData.groupIds(refInd(i)); + if g ~= 0 + members = sort(prepData.refModel.rxns(prepData.groupIds == g)); + rankKeys{i} = members{1}; + end +end + for initStep = 1:length(INITSteps) disp(['ftINIT: Running step ' num2str(initStep)]) stp = INITSteps{initStep}; @@ -259,7 +273,7 @@ params.MIPGap = 0; params.MIPGapAbs = proveAbsGap; try - [deletedRxnsInINIT1, metProduction,fullMipRes,rxnsTurnedOn1,fluxes1] = ftINITInternalAlg(mm,rxnScores,metData,essentialRxns,5,stp.AllowMetSecr,stp.PosRevOff,params, [], fluxes, verbose, 0.1, resolveTies); + [deletedRxnsInINIT1, metProduction,fullMipRes,rxnsTurnedOn1,fluxes1] = ftINITInternalAlg(mm,rxnScores,metData,essentialRxns,5,stp.AllowMetSecr,stp.PosRevOff,params, [], fluxes, verbose, 0.1, resolveTies, rankKeys); fluxes1(toRev) = -fluxes1(toRev); success = true; if fullMipRes.stat ~= 1 @@ -313,7 +327,7 @@ if ~isempty(fullMipRes) startVals = fullMipRes.full; end - [deletedRxnsInINIT1, metProduction,fullMipRes,rxnsTurnedOn1,fluxes1] = ftINITInternalAlg(mm,rxnScores,metData,essentialRxns,5,stp.AllowMetSecr,stp.PosRevOff,params, startVals, fluxes, verbose, 0.1, resolveTies); + [deletedRxnsInINIT1, metProduction,fullMipRes,rxnsTurnedOn1,fluxes1] = ftINITInternalAlg(mm,rxnScores,metData,essentialRxns,5,stp.AllowMetSecr,stp.PosRevOff,params, startVals, fluxes, verbose, 0.1, resolveTies, rankKeys); %This is a bit tricky - since we reversed some reactions, those fluxes also need to be reversed fluxes1(toRev) = -fluxes1(toRev); diff --git a/src/INIT/ftINITInternalAlg.m b/src/INIT/ftINITInternalAlg.m index 47fe78fe..d181505c 100644 --- a/src/INIT/ftINITInternalAlg.m +++ b/src/INIT/ftINITInternalAlg.m @@ -1,4 +1,4 @@ -function [deletedRxns,metProduction,res,turnedOnRxns,fluxes]=ftINITInternalAlg(model,rxnScores,metData,essentialRxns,prodWeight,allowExcretion,remPosRev,params,startVals,fluxes,verbose,forceOn,resolveTies) +function [deletedRxns,metProduction,res,turnedOnRxns,fluxes]=ftINITInternalAlg(model,rxnScores,metData,essentialRxns,prodWeight,allowExcretion,remPosRev,params,startVals,fluxes,verbose,forceOn,resolveTies,rankKeys) % ftINITInternalAlg % This function runs the MILP for a step in ftINIT. % @@ -45,11 +45,17 @@ % free binary, on/off variable (pinned near 1 by the score % objective), so only the negative-score ones are a genuine % tie to resolve. The fluxes are then taken from the -% least-total-flux solution for that on/off state, so the -% directions and forced magnitudes the next ftINIT step -% reads from them are pinned too -- matches raven-toolbox's -% resolve_ties (raven-gecko-parity#104) (optional, default -% false) +% least-total-flux solution for that on/off state and, +% among those, from the one with the least id-rank-weighted +% backward flux, so the directions and forced magnitudes +% the next ftINIT step reads from them are pinned too -- +% matches raven-toolbox's resolve_ties +% (raven-gecko-parity#104) (optional, default false) +% rankKeys cell array with, for each reaction in model.rxns, the key +% that the id-rank phases of resolveTies sort it by +% (optional, default model.rxns). ftINIT passes the smallest +% original reaction id of each merge group, so the order +% does not depend on which member names the merged reaction % % deletedRxns reactions which were deleted by the algorithm (only % rxns included in the problem) @@ -73,6 +79,9 @@ if nargin < 12 || isempty(forceOn) forceOn = 0.1; end +if nargin < 14 || isempty(rankKeys) + rankKeys = model.rxns; +end if nargin < 13 || isempty(resolveTies) resolveTies = false; end @@ -511,7 +520,7 @@ tieProb2b.b = [tieProb.b; kmin+0.5]; tieProb2b.csense = [tieProb.csense 'L']; - [~, sortOrder] = sort(model.rxns(negRxnIdx)); + [~, sortOrder] = sort(rankKeys(negRxnIdx)); ranks = zeros(numel(negRxnIdx),1); ranks(sortOrder) = (1:numel(negRxnIdx))'; idObj = zeros(size(prob.c)); @@ -578,6 +587,43 @@ if isFeasible2c %Only the fluxes change; the status stays that of phase 2b. res.full = res2c.full(1:nCols); + + %Phase 2d: the least total flux can leave a direction open + %when two directions carry the same total flux, for example + %a pair of reversible reactions closing a loop either way. + %Keep the on/off state of phase 2c, hold the total absolute + %flux at most at its phase-2c value plus + %max(1e-7*total, 1e-6), and minimise the backward flux + %max(0,-v) of the reactions whose flux can take either sign, + %weighted by id rank from the end: sorted by rankKeys, the + %first of n reactions gets weight n and the last 1, so a + %reaction with a lower id prefers its forward direction. + if nT > 0 + total = minProb.c' * res2c.full; + dirProb = minProb; + nCols2c = numel(minProb.c); + [~, ord] = sort(rankKeys(tRxns)); + w = zeros(nT,1); + w(ord) = (nT:-1:1)'; + dirProb.A = [minProb.A sparse(size(minProb.A,1),nT); ... + minProb.c' sparse(1,nT); ... + selV sparse(nT,nT) speye(nT)]; + dirProb.a = dirProb.A; + dirProb.b = [minProb.b; total + max(abs(total)*1e-7, 1e-6); zeros(nT,1)]; + dirProb.csense = [minProb.csense 'L' repmat('G',1,nT)]; + dirProb.lb = [minProb.lb; zeros(nT,1)]; + dirProb.ub = [minProb.ub; inf(nT,1)]; + dirProb.vartype = [minProb.vartype repmat('C',1,nT)]; + dirProb.c = [zeros(nCols2c,1); w]; + res2d = optimizeProb(dirProb,minParams,verbose); + [isFeasible2d, isOptimal2d] = checkSolution(res2d); + if ~isOptimal2d + unproven{end+1} = 'phase2d-direction'; + end + if isFeasible2d + res.full = res2d.full(1:nCols); + end + end end end end diff --git a/tests/function_tests/tINIT.m b/tests/function_tests/tINIT.m index cdbd1437..a381bbc4 100644 --- a/tests/function_tests/tINIT.m +++ b/tests/function_tests/tINIT.m @@ -415,6 +415,36 @@ function ftINITInternalAlgResolveTiesGivesMinimalFluxes(testCase) testCase.verifyEqual(fluxes, expected, 'AbsTol', 1e-6); end + function ftINITInternalAlgResolveTiesFixesEqualFluxDirections(testCase) + % The essential E consumes y, which the reversible R1 (y <=> x) + % or R2 (y <=> z) can supply when run backward from the imports + % IN_x or IN_z. Both routes carry the same total flux; with + % resolveTies the higher-ranked reaction carries the backward + % flux: R2 when ranked by id, R1 when rankKeys rank it last. N is + % a negative-score dead end. + testCase.assumeMILPSolver(); + m.rxns = {'IN_x';'IN_z';'R1';'R2';'E';'OUT_p';'N'}; + m.mets = {'x';'y';'z';'p';'d'}; + % INx INz R1 R2 E OUTp N + m.S = sparse([1 0 1 0 0 0 -1; % x + 0 0 -1 -1 -1 0 0; % y + 0 1 0 1 0 0 0; % z + 0 0 0 0 1 -1 0; % p + 0 0 0 0 0 0 1]); % d + m.lb = [0;0;-1000;-1000;0;0;0]; + m.ub = 1000*ones(7,1); + m.rev = double(m.lb < 0); + scores = [0;0;0;0;0;0;-1]; + fluxIn = ones(7,1); + [~,~,~,~,fluxes] = ftINITInternalAlg(m,scores,[],{'E'}, ... + [],false,false,struct(),[],fluxIn,false,0.1,true); + testCase.verifyEqual(fluxes(3:4), [0;-0.1], 'AbsTol', 1e-6); + keys = {'IN_x';'IN_z';'R9';'R2';'E';'OUT_p';'N'}; + [~,~,~,~,fluxes] = ftINITInternalAlg(m,scores,[],{'E'}, ... + [],false,false,struct(),[],fluxIn,false,0.1,true,keys); + testCase.verifyEqual(fluxes(3:4), [-0.1;0], 'AbsTol', 1e-6); + end + function ftINITWithTaskRuns(testCase) % A task requiring e[s] from a[s] forces R2 and R7 to be essential. testCase.assumeMILPSolver(); From ec1a6c196a303524816115a2f9036688fe3d98c9 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 28 Oct 2026 21:22:21 +0100 Subject: [PATCH 04/27] ftINIT: the resolveTies direction phase starts from the minimal-flux solution --- src/INIT/ftINITInternalAlg.m | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/INIT/ftINITInternalAlg.m b/src/INIT/ftINITInternalAlg.m index d181505c..7e8c35f8 100644 --- a/src/INIT/ftINITInternalAlg.m +++ b/src/INIT/ftINITInternalAlg.m @@ -615,6 +615,10 @@ dirProb.ub = [minProb.ub; inf(nT,1)]; dirProb.vartype = [minProb.vartype repmat('C',1,nT)]; dirProb.c = [zeros(nCols2c,1); w]; + %The phase-2c solution is the MIP start: the tight cap on + %the total flux leaves the solver's heuristics little room + %to find a first solution. + dirProb.start = [res2c.full; max(0, -res2c.full(tRxns))]; res2d = optimizeProb(dirProb,minParams,verbose); [isFeasible2d, isOptimal2d] = checkSolution(res2d); if ~isOptimal2d From 33ce66e75acb54ea31c380f0809bf24114162e82 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 28 Oct 2026 23:11:27 +0100 Subject: [PATCH 05/27] mergeCompartments: keep the id of the copy in s or the one that sorts first The merged metabolite was the first copy in model order, so its id depended on the metabolite order of the model. It is now the copy already in compartment s, if any, otherwise the copy whose id sorts first, the rule raven-toolbox's merge_compartments applies. --- src/manipulation/mergeCompartments.m | 31 +++++++++++++++++-------- tests/function_tests/tManipulation.m | 34 ++++++++++++++++++++++++++++ 2 files changed, 55 insertions(+), 10 deletions(-) diff --git a/src/manipulation/mergeCompartments.m b/src/manipulation/mergeCompartments.m index 80c7982e..ce4e80e9 100755 --- a/src/manipulation/mergeCompartments.m +++ b/src/manipulation/mergeCompartments.m @@ -42,8 +42,13 @@ % be useful for example to ensure that there are metabolic capabilities to % synthesize all metabolites. % -% If the metabolite IDs reflect the compartment that they are in the IDs may -% no longer be representative. +% Metabolites with the same name are merged into one. The merged metabolite +% keeps the id and the other fields of one of the copies: the copy that is +% already in compartment "s" if there is one, otherwise the copy whose id +% comes first in character code order (for example ACCOA_c before ACCOA_m). +% With keepUnconstrained only the copies that are not unconstrained are +% candidates. If the metabolite IDs reflect the compartment that they are in +% the IDs may no longer be representative. p=parseRAVENargs(varargin, {'keepUnconstrained',false; 'deleteRxnsWithOneMet',false; 'distReverse',true}); keepUnconstrained=p.keepUnconstrained; @@ -66,9 +71,8 @@ end end -%Loop through each metabolite, and if it is not unconstrained then change -%the S matrix to use the metabolite with the lowest index in model.comps -%instead +%Loop through each metabolite name, and if a copy is not unconstrained then +%change the S matrix to use the copy chosen as described in the Notes %Group metabolite indices by name once, rather than re-scanning all of %model.metNames (via ismember) for every unique name, which is quadratic %in a genome-scale model's (unique names x mets) count. @@ -79,19 +83,26 @@ I=false(numel(model.metNames),1); I(sameNameIdx{i})=true; - %Find the first of those that is not unconstrained. This is the one - %that the other "un-unconstrained" should be changed to. + %Choose the copy that the other "un-unconstrained" copies are merged + %into: one already in compartment "s" if there is one, otherwise the + %one whose id sorts first. if keepUnconstrained==true - mergeTo=find(I & model.unconstrained==false,1); + candidates=find(I & model.unconstrained==false); %This could happen if there is only one metabolite and it is %unconstrained - if isempty(mergeTo) + if isempty(candidates) continue; end else - mergeTo=find(I,1); + candidates=find(I); end + inMerged=candidates(strcmp(model.comps(model.metComps(candidates)),'s')); + if ~isempty(inMerged) + candidates=inMerged; + end + [~,first]=sort(model.mets(candidates)); + mergeTo=candidates(first(1)); I(mergeTo)=false; %Do not do anything for itself I=find(I); diff --git a/tests/function_tests/tManipulation.m b/tests/function_tests/tManipulation.m index e8cfa3b7..aa9ea83f 100644 --- a/tests/function_tests/tManipulation.m +++ b/tests/function_tests/tManipulation.m @@ -510,6 +510,23 @@ function mergeCompartmentsDropsStaleCompMiriams(testCase) end end + function mergeCompartmentsKeepsFirstSortedId(testCase) + % The merged metabolite keeps the id that sorts first, not the + % id of the copy that comes first in the model. + m = tManipulation.mergeToyModel({'ACCOA_m';'ACCOA_c';'X_c'}, {'c';'m'}, [2;1;1]); + evalc('m2 = mergeCompartments(m);'); + testCase.verifyEqual(sort(m2.mets), {'ACCOA_c';'X_c'}); + testCase.verifyEqual(m2.comps, {'s'}); + end + + function mergeCompartmentsPrefersCopyInMergedComp(testCase) + % A copy that is already in compartment s keeps its id, even when + % another copy has an id that sorts before it. + m = tManipulation.mergeToyModel({'ACCOA_c';'ACCOA_z';'X_c'}, {'c';'s'}, [1;2;1]); + evalc('m2 = mergeCompartments(m);'); + testCase.verifyEqual(sort(m2.mets), {'ACCOA_z';'X_c'}); + end + function mergeModelsReturnsStruct(testCase) modelB = removeReactions(testCase.model, testCase.model.rxns(1:5), true, true, true); evalc('merged = mergeModels({testCase.model; modelB});'); @@ -986,6 +1003,23 @@ function expandModelLeavesDnfAlone(testCase) methods (Static, Access = private) + function m = mergeToyModel(mets, comps, metComps) + % Two copies of acetyl-CoA (mets 1 and 2) and X (met 3): + % R1 X => copy 1, R2 copy 2 => X. + m = struct(); + m.id = 'merge'; + m.rxns = {'R1';'R2'}; + m.rxnNames = m.rxns; + m.mets = mets; + m.metNames = {'acetyl-CoA';'acetyl-CoA';'X'}; + m.metComps = metComps; + m.comps = comps; + m.compNames = comps; + m.S = sparse([1 0; 0 -1; -1 1]); + m.lb = [0;0]; m.ub = [1000;1000]; m.rev = [0;0]; m.c = [0;0]; m.b = zeros(3,1); + m.genes = {}; m.grRules = {'';''}; m.rxnGeneMat = sparse(2,0); + end + 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 From 8e38860ee7b5478ac57147029b9af2ba10d0ce70 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 28 Oct 2026 23:11:27 +0100 Subject: [PATCH 06/27] mergeLinear: number groups by their first reaction Group numbers followed the order in which the reactions were merged. They are now 1, 2, ... in the order of the first reaction of each group in origRxnIds, the numbering raven-toolbox's merge_linear uses. --- src/INIT/mergeLinear.m | 10 ++++++++-- tests/function_tests/tINIT.m | 33 +++++++++++++++++++++++++++++++++ 2 files changed, 41 insertions(+), 2 deletions(-) diff --git a/src/INIT/mergeLinear.m b/src/INIT/mergeLinear.m index 068fab65..03d25c4d 100644 --- a/src/INIT/mergeLinear.m +++ b/src/INIT/mergeLinear.m @@ -13,7 +13,8 @@ % reducedModel an updated model structure % origRxnIds Cell array: original rxn ids, used together with the groupIds variable % groupIds Vector of ids: merged rxns have the same id, same length as origRxnIds -% 0 means a reaction was not merged. +% 0 means a reaction was not merged. The groups are numbered 1, 2, ... +% in the order of their first reaction in origRxnIds. % reversedRxns Vector of booleans, saying for each rxn if it is reversed in the combined rxns. % % Usage: [reducedModel,origRxnIds,groupIds]=mergeLinear(model,noMergeRxns) @@ -196,7 +197,12 @@ notInUse=sum(reducedModel.S~=0,2)==0; reducedModel=removeMets(reducedModel,notInUse); end - + +%Number the groups in the order of their first reaction in origRxnIds, so +%the numbers do not depend on the order in which the reactions were merged +merged=groupIds~=0; +[~,~,groupIds(merged)]=unique(groupIds(merged),'stable'); + function [model1,reversedRxns1] = flipRxn(model1, rxnInd, groupIds1, reversedRxns1) model1.S(:,rxnInd) = model1.S(:,rxnInd)*-1; %swap the bounds diff --git a/tests/function_tests/tINIT.m b/tests/function_tests/tINIT.m index a381bbc4..e5dbcc46 100644 --- a/tests/function_tests/tINIT.m +++ b/tests/function_tests/tINIT.m @@ -179,6 +179,39 @@ function mergeLinearGroupsAndScores(testCase) testCase.verifyTrue(all(find(reversedRxns) == [6;9])); end + function mergeLinearNumbersGroupsByFirstRxn(testCase) + % Group numbers follow the first reaction of each group in + % model.rxns, not the merge order. The metabolites are visited + % in model.mets order, so c (R5, R6) merges first. + model = struct(); + model.id = 'numbering'; + model.rxns = {'R1';'R2';'R3';'R4';'R5';'R6'}; + model.rxnNames = model.rxns; + model.mets = {'z';'c';'w';'x';'a';'m';'b';'y'}; + model.metNames = model.mets; + model.comps = {'c'}; + model.compNames = {'c'}; + model.metComps = ones(8,1); + subs = {'x';'a';'m';'b';'z';'c'}; + prods = {'a';'m';'b';'y';'c';'w'}; + model.S = sparse(8,6); + for i = 1:6 + model.S(strcmp(model.mets,subs{i}),i) = -1; + model.S(strcmp(model.mets,prods{i}),i) = 1; + end + model.lb = zeros(6,1); + model.ub = 1000*ones(6,1); + model.rev = zeros(6,1); + model.c = zeros(6,1); + model.b = zeros(8,1); + model.genes = cell(0,1); + model.grRules = repmat({''},6,1); + model.rxnGeneMat = sparse(6,0); + evalc('[reducedModel,~,groupIds] = mergeLinear(model, {});'); + testCase.verifyEqual(groupIds(:), [1;1;1;1;2;2]); + testCase.verifyEqual(sort(reducedModel.rxns), {'R1';'R5'}); + end + function scoreModelRuns(testCase) arrayData.genes = testCase.model.genes; arrayData.tissues = {'t1'}; From 629d56be10a872dd11a2a9069b4083edd90cbc44 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 28 Oct 2026 23:12:58 +0100 Subject: [PATCH 07/27] grRuleToDNF: help states the parseGrRule grammar and the RAVEN:badGrRule error; tUtils tests operator spellings and malformed rules --- src/utils/grRuleToDNF.m | 9 ++++++++- tests/function_tests/tUtils.m | 15 +++++++++++++++ 2 files changed, 23 insertions(+), 1 deletion(-) diff --git a/src/utils/grRuleToDNF.m b/src/utils/grRuleToDNF.m index 34da5a64..f2ea4465 100644 --- a/src/utils/grRuleToDNF.m +++ b/src/utils/grRuleToDNF.m @@ -8,7 +8,10 @@ % Parameters % ---------- % rule : char or string or struct -% a grRule, or a tree from parseGrRule. +% a grRule, or a tree from parseGrRule. A grRule is read by parseGrRule: +% genes joined by "and"/"or" in any case or by "&", "&&", "|", "||", +% each operator separated from the gene ids by whitespace or a bracket, +% with AND binding tighter than OR. % % Returns % ------- @@ -18,6 +21,10 @@ % % Notes % ----- +% A malformed rule (an operator without an operand, unbalanced or empty +% parentheses, two genes without an operator between them) raises +% RAVEN:badGrRule. +% % Gene order within a clause follows the order the genes appear in the rule. % Duplicates are not removed: "G1 and G1" stays a two-element clause, since % collapsing it would silently rewrite the user's rule. diff --git a/tests/function_tests/tUtils.m b/tests/function_tests/tUtils.m index 858ed0b2..1b45ab53 100644 --- a/tests/function_tests/tUtils.m +++ b/tests/function_tests/tUtils.m @@ -181,6 +181,21 @@ function grRuleToDNFDistributes(testCase) testCase.verifyEqual(numel(grRuleToDNF('(g1 or g2) and (g3 or g4)')), 4); end + function grRuleToDNFTakesTheParseGrRuleGrammar(testCase) + % Doubled symbols and operators in any case are operators; an + % operator glued to a gene id is part of that id. + testCase.verifyEqual(grRuleToDNF('G1 && (G2 || G3)'), {{'G1','G2'},{'G1','G3'}}); + testCase.verifyEqual(grRuleToDNF('G1 Or G2 aNd G3'), {{'G1'},{'G2','G3'}}); + testCase.verifyEqual(grRuleToDNF('G1&G2 or G3'), {{'G1&G2'},{'G3'}}); + end + + function grRuleToDNFRejectsMalformed(testCase) + rules = {'G1 and', '(G1 or G2', 'G1 G2', '()', 'G1 or or G2', 'G1 xor G2'}; + for k = 1:numel(rules) + testCase.verifyError(@() grRuleToDNF(rules{k}), 'RAVEN:badGrRule', rules{k}); + end + end + function grRuleToDNFSimpleCases(testCase) testCase.verifyEqual(grRuleToDNF(''), {}); testCase.verifyEqual(grRuleToDNF('G1'), {{'G1'}}); From 62717488c839c336768a7f0903c1e8920ac22dca Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 28 Oct 2026 23:12:58 +0100 Subject: [PATCH 08/27] gapFillTopological: help states that the default targets are every metabolite of the objective reaction; tGapfilling tests the default targets --- src/gapfilling/gapFillTopological.m | 7 ++++--- tests/function_tests/tGapfilling.m | 15 +++++++++++++++ 2 files changed, 19 insertions(+), 3 deletions(-) diff --git a/src/gapfilling/gapFillTopological.m b/src/gapfilling/gapFillTopological.m index fe6e335e..e5c36a2a 100644 --- a/src/gapfilling/gapFillTopological.m +++ b/src/gapfilling/gapFillTopological.m @@ -34,8 +34,9 @@ % metabolites of every exchange reaction whose bounds let it supply the % metabolite, i.e. getExchangeRxns 'uptake' plus 'reverse'. % 'targets' (default []) -% Metabolite IDs (cell array) that should be produced. Default: -% substrates of the objective (biomass) reaction. +% Metabolite IDs (cell array) that should be produced. Default: every +% metabolite of the objective (biomass) reaction, substrates and +% products; the first one when there are several objective reactions. % 'verbose' (default true) % print a summary of the results. % @@ -99,7 +100,7 @@ % ---- Identify target metabolites ---- if isempty(targets) - % Default: substrates of the objective (biomass) reaction + % Default: every metabolite of the objective (biomass) reaction objRxns = find(model.c ~= 0); if isempty(objRxns) error(['gapFillTopological: no objective function defined. ' ... diff --git a/tests/function_tests/tGapfilling.m b/tests/function_tests/tGapfilling.m index acd68455..9a4c2ee0 100644 --- a/tests/function_tests/tGapfilling.m +++ b/tests/function_tests/tGapfilling.m @@ -350,6 +350,21 @@ function gapFillTopologicalBlocksWithoutExchange(testCase) testCase.verifyNumElements(result.blockedMets, 3); end + function gapFillTopologicalDefaultTargetsAreObjectiveMets(testCase) + % Objective R2 (B -> C): the default targets are its substrate B + % and its product C. Without an exchange both are blocked. + m = testCase.chainModel(); + evalc('result = gapFillTopological(m, m, ''verbose'', false);'); + testCase.verifyEqual(sort(result.blockedMets(:)), {'B';'C'}); + end + + function gapFillTopologicalDefaultTargetsNeedObjective(testCase) + m = testCase.chainModel(); + m.c = [0;0]; + testCase.verifyError(@() gapFillTopological(m, m, 'seeds', {'A'}, ... + 'verbose', false), ?MException); + end + function gapFillMILPRepairsGrowth(testCase) testCase.assumeMILPSolver(); modelDB = testCase.model; modelDB.id = 'DB'; From 7e2000533e6daf3cae3cf4029bde38644b980a3f Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 28 Oct 2026 23:52:01 +0100 Subject: [PATCH 09/27] gapFillTopological: the no-objective error is RAVEN:badInput Without targets and without an objective reaction the function raised an error without identifier, which tGapfilling could only test as any MException. It now raises RAVEN:badInput, which the help states and tGapfilling tests. --- src/gapfilling/gapFillTopological.m | 4 +++- tests/function_tests/tGapfilling.m | 2 +- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/src/gapfilling/gapFillTopological.m b/src/gapfilling/gapFillTopological.m index e5c36a2a..bf832c02 100644 --- a/src/gapfilling/gapFillTopological.m +++ b/src/gapfilling/gapFillTopological.m @@ -37,6 +37,8 @@ % Metabolite IDs (cell array) that should be produced. Default: every % metabolite of the objective (biomass) reaction, substrates and % products; the first one when there are several objective reactions. +% Without targets and without an objective reaction the function raises +% RAVEN:badInput. % 'verbose' (default true) % print a summary of the results. % @@ -103,7 +105,7 @@ % Default: every metabolite of the objective (biomass) reaction objRxns = find(model.c ~= 0); if isempty(objRxns) - error(['gapFillTopological: no objective function defined. ' ... + error('RAVEN:badInput', ['gapFillTopological: no objective function defined. ' ... 'Provide targets manually via the ''targets'' option.']); end biomassIdx = objRxns(1); diff --git a/tests/function_tests/tGapfilling.m b/tests/function_tests/tGapfilling.m index 9a4c2ee0..b4492fa6 100644 --- a/tests/function_tests/tGapfilling.m +++ b/tests/function_tests/tGapfilling.m @@ -362,7 +362,7 @@ function gapFillTopologicalDefaultTargetsNeedObjective(testCase) m = testCase.chainModel(); m.c = [0;0]; testCase.verifyError(@() gapFillTopological(m, m, 'seeds', {'A'}, ... - 'verbose', false), ?MException); + 'verbose', false), 'RAVEN:badInput'); end function gapFillMILPRepairsGrowth(testCase) From be2d17b8cb3c4de1fa87006eaabeec603a789f5d Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 28 Oct 2026 23:52:01 +0100 Subject: [PATCH 10/27] tUtils: grRuleToDNF splits words at the characters isspace counts char(28) and char(8195) separate words; char(133), char(160), char(8199) and char(8239), for which isspace is false, are part of a word. --- tests/function_tests/tUtils.m | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/tests/function_tests/tUtils.m b/tests/function_tests/tUtils.m index 1b45ab53..b85024ca 100644 --- a/tests/function_tests/tUtils.m +++ b/tests/function_tests/tUtils.m @@ -189,6 +189,17 @@ function grRuleToDNFTakesTheParseGrRuleGrammar(testCase) testCase.verifyEqual(grRuleToDNF('G1&G2 or G3'), {{'G1&G2'},{'G3'}}); end + function grRuleToDNFWhitespaceIsIsspace(testCase) + % isspace is true for char(28) and char(8195) and false for + % the no-break space char(160) and char(8239), which are then + % part of a word. + testCase.verifyEqual(grRuleToDNF(['G1' char(28) 'or' char(8195) 'G2']), {{'G1'},{'G2'}}); + for c = [133 160 8199 8239] + testCase.verifyError(@() grRuleToDNF(['G1' char(c) 'and G2']), 'RAVEN:badGrRule'); + testCase.verifyEqual(grRuleToDNF(['G1' char(c) 'G2']), {{['G1' char(c) 'G2']}}); + end + end + function grRuleToDNFRejectsMalformed(testCase) rules = {'G1 and', '(G1 or G2', 'G1 G2', '()', 'G1 or or G2', 'G1 xor G2'}; for k = 1:numel(rules) From 52f807f940da29e671e0add19e73adabe8a34aa5 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 28 Oct 2026 23:52:01 +0100 Subject: [PATCH 11/27] CHANGELOG: Unreleased section for mergeCompartments and mergeLinear --- CHANGELOG.md | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 2e5acd85..0c65bc73 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,14 @@ ## Unreleased +### Manipulation + +* `mergeCompartments`: a merged metabolite is the copy already in compartment `s`, if any, + otherwise the copy whose id sorts first, with every field of that copy (was: the first copy + in model order), so `ACCOA_c` and `ACCOA_m` become `ACCOA_c` whatever their order. + +### Context-specific models (ftINIT / tINIT) + * `ftINIT` with `resolveTies`: each main-extraction step takes its fluxes from the least-total-flux solution for the tie-broken on/off state, so the direction and magnitude at which the next step forces a reaction no longer depend on the solver's vertex. @@ -11,6 +19,8 @@ direction of a flux that carries the same total either way no longer depends on the solver. The id-rank phases rank a merged reaction by the smallest original reaction id in its group, so they no longer depend on which member names it. +* `mergeLinear`: groups are numbered 1, 2, ... in the order of the first reaction of each + group in `origRxnIds` (was: the order in which the reactions were merged). ## 3.0.0b1 — 2026-09-07 From 23bd87be5ec739c5a61d5835f970822a27e42409 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 28 Oct 2026 23:52:01 +0100 Subject: [PATCH 12/27] CHANGELOG: Unreleased section for grRuleToDNF and gapFillTopological --- CHANGELOG.md | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 0c65bc73..c5c4755e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,12 @@ otherwise the copy whose id sorts first, with every field of that copy (was: the first copy in model order), so `ACCOA_c` and `ACCOA_m` become `ACCOA_c` whatever their order. +### Gap-filling + +* `gapFillTopological`: without `targets` and without an objective reaction it raises + `RAVEN:badInput` (was: an error without identifier). Help states the default seeds and + targets. + ### Context-specific models (ftINIT / tINIT) * `ftINIT` with `resolveTies`: each main-extraction step takes its fluxes from the @@ -22,6 +28,12 @@ * `mergeLinear`: groups are numbered 1, 2, ... in the order of the first reaction of each group in `origRxnIds` (was: the order in which the reactions were merged). +### Utilities + +* `grRuleToDNF`: help states the grammar it reads, that of `parseGrRule` (`and`/`or` in any + case, `&`, `&&`, `|`, `||`, operators as whole words, so `G1&G2` is one gene id), and that a + malformed rule raises `RAVEN:badGrRule`. + ## 3.0.0b1 — 2026-09-07 First beta toward 3.0. Changes below are relative to RAVEN 2.9.3, the last 2.x release: From 66c2df640b2d147537ed5efa1fd74e4da25f4003 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 29 Oct 2026 00:26:03 +0100 Subject: [PATCH 13/27] mergeLinear: a reaction flipped in a later pass is marked reversed by its id flipRxn indexed groupIds and reversedRxns, which follow origRxnIds, with the index of the reduced model, so after the first pass it marked the reaction at the same original index instead of the flipped one. --- CHANGELOG.md | 2 ++ src/INIT/mergeLinear.m | 11 +++++++---- tests/function_tests/tINIT.m | 36 ++++++++++++++++++++++++++++++++++++ 3 files changed, 45 insertions(+), 4 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index c5c4755e..c0374ffd 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -27,6 +27,8 @@ its group, so they no longer depend on which member names it. * `mergeLinear`: groups are numbered 1, 2, ... in the order of the first reaction of each group in `origRxnIds` (was: the order in which the reactions were merged). +* `mergeLinear`: a reaction flipped after the first merge pass is the one marked in + `reversedRxns` (was: the reaction at the same index of the original model). ### Utilities diff --git a/src/INIT/mergeLinear.m b/src/INIT/mergeLinear.m index 03d25c4d..66847e98 100644 --- a/src/INIT/mergeLinear.m +++ b/src/INIT/mergeLinear.m @@ -214,10 +214,13 @@ %now take care of the reversedRxns - if this is a group, reverse all of the %reactions in the group in the reversedRxns index - they will all be reversed at the - %same time since they are the same rxn. - rxnIndices = rxnInd; - if groupIds1(rxnInd) > 0 - rxnIndices = find(groupIds1 == groupIds1(rxnInd)); + %same time since they are the same rxn. rxnInd indexes model1, which loses + %reactions after each pass, while groupIds1 and reversedRxns1 follow + %origRxnIds, so the reaction is looked up by its id. + origInd = origRxnIdMap(model1.rxns{rxnInd}); + rxnIndices = origInd; + if groupIds1(origInd) > 0 + rxnIndices = find(groupIds1 == groupIds1(origInd)); end reversedRxns1(rxnIndices) = ~reversedRxns1(rxnIndices); end diff --git a/tests/function_tests/tINIT.m b/tests/function_tests/tINIT.m index e5dbcc46..fe4fe7dd 100644 --- a/tests/function_tests/tINIT.m +++ b/tests/function_tests/tINIT.m @@ -212,6 +212,42 @@ function mergeLinearNumbersGroupsByFirstRxn(testCase) testCase.verifyEqual(sort(reducedModel.rxns), {'R1';'R5'}); end + function mergeLinearFlipsRxnInLaterPass(testCase) + % R1 and R2 merge through A in the first pass and R2 is removed, + % so R4 sits at index 3 of the reduced model. The second pass + % merges R1 and R4 through Y and flips R4, which must be the + % reaction marked as reversed (not R3, the third original one). + model = struct(); + model.id = 'lateflip'; + model.rxns = {'R1';'R2';'R3';'R4';'R5';'R6'}; + model.rxnNames = model.rxns; + model.mets = {'A';'Y';'Q';'W'}; + model.metNames = model.mets; + model.comps = {'c'}; + model.compNames = {'c'}; + model.metComps = ones(4,1); + S = zeros(4,6); + S([1 2],1) = [1;1]; % R1: -> A + Y + S([1 2],2) = [-1;1]; % R2: A -> Y + S(3,3) = -1; % R3: Q -> + S([4 2],4) = [-1;1]; % R4: W <=> Y + S(4,5) = 1; % R5: -> W + S(4,6) = -1; % R6: W -> + model.S = sparse(S); + model.lb = [0;0;0;-1000;0;0]; + model.ub = 1000*ones(6,1); + model.rev = double(model.lb < 0); + model.c = zeros(6,1); + model.b = zeros(4,1); + model.genes = cell(0,1); + model.grRules = repmat({''},6,1); + model.rxnGeneMat = sparse(6,0); + evalc('[reducedModel,origRxnIds,groupIds,reversedRxns] = mergeLinear(model, {});'); + testCase.verifyEqual(groupIds(:), [1;1;0;1;0;0]); + testCase.verifyEqual(origRxnIds(reversedRxns), {'R4'}); + testCase.verifyEqual(sort(reducedModel.rxns), {'R1';'R3';'R5';'R6'}); + end + function scoreModelRuns(testCase) arrayData.genes = testCase.model.genes; arrayData.tissues = {'t1'}; From c8daa8527cc960c627e57006d5f5272cf010f366 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 29 Oct 2026 00:26:03 +0100 Subject: [PATCH 14/27] removeGenes: help states that genes no remaining reaction uses are removed --- CHANGELOG.md | 1 + src/manipulation/removeGenes.m | 6 +++++- 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index c0374ffd..87599882 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,7 @@ * `mergeCompartments`: a merged metabolite is the copy already in compartment `s`, if any, otherwise the copy whose id sorts first, with every field of that copy (was: the first copy in model order), so `ACCOA_c` and `ACCOA_m` become `ACCOA_c` whatever their order. +* `removeGenes`: the help states that genes no remaining reaction uses are removed as well. ### Gap-filling diff --git a/src/manipulation/removeGenes.m b/src/manipulation/removeGenes.m index e2a50f32..c6629859 100755 --- a/src/manipulation/removeGenes.m +++ b/src/manipulation/removeGenes.m @@ -23,7 +23,11 @@ % Returns % ------- % reducedModel : struct -% an updated model structure. +% an updated model structure. Besides genesToRemove, every gene that no +% remaining reaction used before the removal is deleted as well, also +% when it was not requested (for example a gene that no reaction uses). +% A gene whose only reaction is kept but blocked (removeBlockedRxns +% false), such as the other subunit of a complex that lost one, stays. % % Examples % -------- From 95e47b2cc171ac2737f7dab6134a789bac9f1c0f Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 29 Oct 2026 01:12:09 +0100 Subject: [PATCH 15/27] removeGenes: help states that nothing changes when no requested gene is in the model --- CHANGELOG.md | 4 +++- src/manipulation/removeGenes.m | 14 +++++++++----- 2 files changed, 12 insertions(+), 6 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 87599882..ea731448 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,7 +7,9 @@ * `mergeCompartments`: a merged metabolite is the copy already in compartment `s`, if any, otherwise the copy whose id sorts first, with every field of that copy (was: the first copy in model order), so `ACCOA_c` and `ACCOA_m` become `ACCOA_c` whatever their order. -* `removeGenes`: the help states that genes no remaining reaction uses are removed as well. +* `removeGenes`: the help states that, when at least one of the requested genes is in the model, + genes no remaining reaction uses are removed as well; when none is, no gene, reaction or bound + changes. ### Gap-filling diff --git a/src/manipulation/removeGenes.m b/src/manipulation/removeGenes.m index c6629859..00f3e30f 100755 --- a/src/manipulation/removeGenes.m +++ b/src/manipulation/removeGenes.m @@ -23,11 +23,15 @@ % Returns % ------- % reducedModel : struct -% an updated model structure. Besides genesToRemove, every gene that no -% remaining reaction used before the removal is deleted as well, also -% when it was not requested (for example a gene that no reaction uses). -% A gene whose only reaction is kept but blocked (removeBlockedRxns -% false), such as the other subunit of a complex that lost one, stays. +% an updated model structure. Genes in genesToRemove that are not in the +% model are ignored. When none of them is in the model, no gene, reaction +% or bound changes (with standardizeRules true the gene rules are still +% written in the standard format). Otherwise, besides genesToRemove, every +% gene that no remaining reaction used before the removal is deleted as +% well, also when it was not requested (for example a gene that no +% reaction uses). A gene whose only reaction is kept but blocked +% (removeBlockedRxns false), such as the other subunit of a complex that +% lost one, stays. % % Examples % -------- From 16b387b4efacd5759be16e7f3dda71557152d077 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 29 Oct 2026 01:12:09 +0100 Subject: [PATCH 16/27] tINIT: ftINIT ranks a merged reaction by the smallest original id of its group --- tests/function_tests/tINIT.m | 64 ++++++++++++++++++++++++++++++++++++ 1 file changed, 64 insertions(+) diff --git a/tests/function_tests/tINIT.m b/tests/function_tests/tINIT.m index fe4fe7dd..2486b126 100644 --- a/tests/function_tests/tINIT.m +++ b/tests/function_tests/tINIT.m @@ -453,6 +453,34 @@ function ftINITResolveTiesAndProveAbsGapStillRunCorrectly(testCase) testCase.verifyTrue(all(strcmp(resBoth.rxns, expected))); end + function ftINITResolveTiesRanksMergedRxnBySmallestOriginalId(testCase) + % a[c] reaches m[c] by path P (Rm, Rc) or path Q (Rn, Rb), each + % path one linear merge group with equal negative scores. Model + % order and group numbers put P first, and the merged reactions + % are named by their first members Rm and Rn, but the smallest + % original id of Q (Rb) sorts before that of P (Rc). resolveTies + % keeps Q; ranking by the merged id, the model order or the group + % number would keep P. + testCase.assumeMILPSolver(); + testModel = getTstModelTwoPaths(); + evalc('prepData = prepINITModel(testModel, getTstModelTasks(), {}, false, {}, ''s'');'); + g = @(r) prepData.groupIds(strcmp(prepData.refModel.rxns, r)); + testCase.verifyEqual(g('Rm'), g('Rc')); + testCase.verifyEqual(g('Rn'), g('Rb')); + testCase.verifyLessThan(g('Rm'), g('Rn')); + testCase.verifyTrue(all(ismember({'Rm';'Rn'}, prepData.minModel.rxns))); + + arrayData.genes = testModel.genes; + arrayData.tissues = {'a'}; + arrayData.levels = getExprForRxnScore([-1;-1;-1;-1;5;5]); + arrayData.threshold = 1; + evalc(['resModel = ftINIT(prepData,arrayData.tissues{1},[],[],' ... + 'arrayData,[],getINITSteps(),true,true,struct(),false,' ... + '''resolveTies'',true);']); + testCase.verifyTrue(all(ismember({'Rn';'Rb'}, resModel.rxns))); + testCase.verifyFalse(any(ismember({'Rm';'Rc'}, resModel.rxns))); + end + function ftINITInternalAlgResolveTiesGivesMinimalFluxes(testCase) % With resolveTies the fluxes come from the least-total-flux % solution for the tie-broken on/off state. The positive @@ -946,6 +974,42 @@ function scoreModelRejectsUnknownDataPrecedence(testCase) testModel.b = repmat(0,8,1); end +function testModel = getTstModelTwoPaths() +%a[c] reaches m[c] by Rm, Rc or by Rn, Rb; E, R7 and R8 export e[s]. +testModel = struct(); +testModel.id = 'twoPaths'; +testModel.rxns = {}; +testModel.S = []; +testModel.rev = []; +testModel.mets = {'as';'ac';'xc';'yc';'mc';'ec';'es'}; +testModel.metNames = {'a';'a';'x';'y';'m';'e';'e'}; +testModel.comps = {'s';'c'}; +testModel.compNames = testModel.comps; +testModel.metComps = [1;2;2;2;2;2;1]; +testModel.genes = {'Gm';'Gc';'Gn';'Gb';'GE';'G7'}; +testModel.grRules = {}; +testModel.rxnGeneMat = []; + +rxnsToAdd = struct(); +rxnsToAdd.rxns = {'R1';'R2';'Rm';'Rc';'Rn';'Rb';'E';'R7';'R8'}; +rxnsToAdd.equations = {'=> a[s]';... + 'a[s] <=> a[c]';... + 'a[c] => x[c]';... + 'x[c] => m[c]';... + 'a[c] => y[c]';... + 'y[c] => m[c]';... + 'm[c] => e[c]';... + 'e[c] => e[s]';... + 'e[s] =>'}; +rxnsToAdd.grRules = {'';'';'Gm';'Gc';'Gn';'Gb';'GE';'G7';''}; +evalc('testModel = addRxns(testModel,rxnsToAdd, 3);'); +testModel.c = [0;0;0;0;0;0;0;0;1]; +testModel.ub = repmat(1000,9,1); +testModel.lb = [0;-1000;0;0;0;0;0;0;0]; +testModel.rxnNames = testModel.rxns; +testModel.b = zeros(7,1); +end + function testModelRxnScores = getTstModelRxnScores() testModelRxnScores = [-2;-2;-1;7;0.5;0.5;-1;-2;-3;3.5]; end From 052263add161a8d15250db586ef51b8e9ea7f9e9 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 29 Oct 2026 01:25:01 +0100 Subject: [PATCH 17/27] tINIT: the paths of the rank-key test model are reversible, as prepINITModel needs --- tests/function_tests/tINIT.m | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/tests/function_tests/tINIT.m b/tests/function_tests/tINIT.m index 2486b126..4f00e13a 100644 --- a/tests/function_tests/tINIT.m +++ b/tests/function_tests/tINIT.m @@ -975,7 +975,9 @@ function scoreModelRejectsUnknownDataPrecedence(testCase) end function testModel = getTstModelTwoPaths() -%a[c] reaches m[c] by Rm, Rc or by Rn, Rb; E, R7 and R8 export e[s]. +%a[c] reaches m[c] by Rm, Rc or by Rn, Rb; E, R7 and R8 export e[s]. The +%paths are reversible: simplifyModel in prepINITModel fails on a model +%without a reversible reaction after the first simplification. testModel = struct(); testModel.id = 'twoPaths'; testModel.rxns = {}; @@ -994,10 +996,10 @@ function scoreModelRejectsUnknownDataPrecedence(testCase) rxnsToAdd.rxns = {'R1';'R2';'Rm';'Rc';'Rn';'Rb';'E';'R7';'R8'}; rxnsToAdd.equations = {'=> a[s]';... 'a[s] <=> a[c]';... - 'a[c] => x[c]';... - 'x[c] => m[c]';... - 'a[c] => y[c]';... - 'y[c] => m[c]';... + 'a[c] <=> x[c]';... + 'x[c] <=> m[c]';... + 'a[c] <=> y[c]';... + 'y[c] <=> m[c]';... 'm[c] => e[c]';... 'e[c] => e[s]';... 'e[s] =>'}; @@ -1005,7 +1007,7 @@ function scoreModelRejectsUnknownDataPrecedence(testCase) evalc('testModel = addRxns(testModel,rxnsToAdd, 3);'); testModel.c = [0;0;0;0;0;0;0;0;1]; testModel.ub = repmat(1000,9,1); -testModel.lb = [0;-1000;0;0;0;0;0;0;0]; +testModel.lb = [0;-1000;-1000;-1000;-1000;-1000;0;0;0]; testModel.rxnNames = testModel.rxns; testModel.b = zeros(7,1); end From bd238a8cc695596532e2ed75e5af4f787720b709 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 29 Oct 2026 07:34:36 +0100 Subject: [PATCH 18/27] ftINITInternalAlg: resolveTies skips the tie-break after a time-limited primary solve and gives phases 2a and 2b half of TimeLimit --- CHANGELOG.md | 4 ++++ src/INIT/ftINIT.m | 6 +++++- src/INIT/ftINITInternalAlg.m | 18 ++++++++++++++++-- 3 files changed, 25 insertions(+), 3 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index ea731448..525e2813 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -28,6 +28,10 @@ direction of a flux that carries the same total either way no longer depends on the solver. The id-rank phases rank a merged reaction by the smallest original reaction id in its group, so they no longer depend on which member names it. +* `ftINITInternalAlg` with `resolveTies`: when the primary solve stops at the time limit, the + tie-break phases are skipped and the primary solution is kept, as in raven-toolbox (was: the + phases ran on the unproven incumbent). The count and id-rank phases get half of + `params.TimeLimit` each, rounded down and at least 1 s (was: the full limit). * `mergeLinear`: groups are numbered 1, 2, ... in the order of the first reaction of each group in `origRxnIds` (was: the order in which the reactions were merged). * `mergeLinear`: a reaction flipped after the first merge pass is the one marked in diff --git a/src/INIT/ftINIT.m b/src/INIT/ftINIT.m index 6efa2e47..f8449204 100644 --- a/src/INIT/ftINIT.m +++ b/src/INIT/ftINIT.m @@ -79,7 +79,11 @@ % flux (a merged reaction ranked by the smallest original id in its % group). Applies to every main-extraction step and to the final % gap-filling step (fitTasks' resolveTies) -- see raven-gecko-parity#104 -% and raven-toolbox's resolve_ties (default false). +% and raven-toolbox's resolve_ties. A step whose primary solve stops at +% the time limit keeps its primary solution without these phases; the +% count and id-rank phases get half of the step's TimeLimit each, and a warning +% names each phase that does not prove optimal or finds no solution +% (default false). % proveAbsGap : double % if set, prove each step to this fixed *absolute* MILP gap in a single % solve, instead of the relative-gap escalation below (stp.MILPParams). diff --git a/src/INIT/ftINITInternalAlg.m b/src/INIT/ftINITInternalAlg.m index 7e8c35f8..6caccf0e 100644 --- a/src/INIT/ftINITInternalAlg.m +++ b/src/INIT/ftINITInternalAlg.m @@ -50,7 +50,14 @@ % backward flux, so the directions and forced magnitudes % the next ftINIT step reads from them are pinned too -- % matches raven-toolbox's resolve_ties -% (raven-gecko-parity#104) (optional, default false) +% (raven-gecko-parity#104). The phases are skipped when +% the primary solve stops at params.TimeLimit, and the +% primary solution is kept. Phases 2a and 2b get half of +% params.TimeLimit each (rounded down, at least 1 s), the +% flux phases the full limit. A warning names each phase +% that does not prove optimal or finds no solution; when a +% flux phase finds none, the fluxes of the phase before it +% are kept (optional, default false) % rankKeys cell array with, for each reaction in model.rxns, the key % that the id-rank phases of resolveTies sort it by % (optional, default model.rxns). ftINIT passes the smallest @@ -465,7 +472,10 @@ error('RAVEN:badInput', '%s', EM); end -if resolveTies && (nNegIrrev + nNegRev) > 0 +if resolveTies && (nNegIrrev + nNegRev) > 0 && ~res.hitTimeLimit + %Skipped when the primary solve stopped at the time limit: the floor + %constraint would pin an unproven incumbent, so the phases cannot select + %among true optima. %This MILP is highly degenerate: many reaction subsets reach the same %score optimum, and the solver returns an arbitrary one -- reproducible %for a fixed solver build + Seed, but fragile to a solver version bump @@ -495,6 +505,10 @@ tieParams = params; tieParams.MIPGap = 0; tieParams.MIPGapAbs = 0.4; + %Phases 2a and 2b get half of the time limit each, as in raven-toolbox. + if isfield(tieParams,'TimeLimit') + tieParams.TimeLimit = max(1,floor(tieParams.TimeLimit*0.5)); + end unproven = {}; From 202f87ad495c45153bcf4dd0e97ce62b41f2d35e Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 29 Oct 2026 07:34:36 +0100 Subject: [PATCH 19/27] tINIT: resolveTies runs no tie-break phase after a time-limited primary solve --- tests/function_tests/tINIT.m | 39 ++++++++++++++++++++++++++++++++++++ 1 file changed, 39 insertions(+) diff --git a/tests/function_tests/tINIT.m b/tests/function_tests/tINIT.m index 4f00e13a..4c4e9c32 100644 --- a/tests/function_tests/tINIT.m +++ b/tests/function_tests/tINIT.m @@ -542,6 +542,45 @@ function ftINITInternalAlgResolveTiesFixesEqualFluxDirections(testCase) testCase.verifyEqual(fluxes(3:4), [-0.1;0], 'AbsTol', 1e-6); end + function ftINITInternalAlgResolveTiesSkippedAfterTimeLimit(testCase) + % When the primary solve stops at the time limit, resolveTies + % keeps the primary solution and runs no tie-break phase. A + % wrapper of optimizeProb earlier on the path counts the solves + % and reports each as stopped at the time limit. + testCase.assumeMILPSolver(); + global TINIT_OPTIMIZEPROB TINIT_OPTIMIZEPROB_CALLS + TINIT_OPTIMIZEPROB = @optimizeProb; + TINIT_OPTIMIZEPROB_CALLS = 0; + d = tempname; mkdir(d); + fid = fopen(fullfile(d,'optimizeProb.m'),'w'); + fprintf(fid,'%s\n', ... + 'function res = optimizeProb(varargin)', ... + 'global TINIT_OPTIMIZEPROB TINIT_OPTIMIZEPROB_CALLS', ... + 'res = TINIT_OPTIMIZEPROB(varargin{:});', ... + 'TINIT_OPTIMIZEPROB_CALLS = TINIT_OPTIMIZEPROB_CALLS + 1;', ... + 'res.hitTimeLimit = true;', ... + 'end'); + fclose(fid); + testCase.addTeardown(@() rmdir(d,'s')); + addpath(d,'-begin'); + testCase.addTeardown(@() rmpath(d)); + testCase.addTeardown(@() clear('global','TINIT_OPTIMIZEPROB','TINIT_OPTIMIZEPROB_CALLS')); + % R1 and R2 can each feed the essential E: a tie for resolveTies. + m.rxns = {'R1';'R2';'E';'EX_a';'EX_p'}; + m.mets = {'a';'m';'p'}; + m.S = sparse([-1 -1 0 -1 0; + 1 1 -1 0 0; + 0 0 1 0 -1]); + m.lb = [0;0;0;-1000;-1000]; + m.ub = 1000*ones(5,1); + m.rev = double(m.lb < 0); + [~,~,res,turnedOn] = ftINITInternalAlg(m,[-1;-1;0;0;0],[],{'E'}, ... + [],false,false,struct(),[],0.1*ones(5,1),false,0.1,true); + testCase.verifyEqual(TINIT_OPTIMIZEPROB_CALLS, 1); + testCase.verifyTrue(res.hitTimeLimit); + testCase.verifyEqual(nnz(turnedOn(1:2)), 1); + end + function ftINITWithTaskRuns(testCase) % A task requiring e[s] from a[s] forces R2 and R7 to be essential. testCase.assumeMILPSolver(); From c3635e807f2882c7f1dc5b1670c6decda6831cff Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 29 Oct 2026 07:43:24 +0100 Subject: [PATCH 20/27] prepINITModel: prepData carries reversedRxns of mergeLinear --- src/INIT/prepINITModel.m | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/INIT/prepINITModel.m b/src/INIT/prepINITModel.m index e129778c..d2a27ed3 100644 --- a/src/INIT/prepINITModel.m +++ b/src/INIT/prepINITModel.m @@ -194,7 +194,7 @@ %Now, merge all linear dependent rxns - this reduces the size of the model a lot. %Size goes from 11888 rxns to 7922 rxns for human-GEM, so this is an important step -[minModel3,origRxnIds,groupIds]=mergeLinear(minModel2, {}); +[minModel3,origRxnIds,groupIds,reversedRxns]=mergeLinear(minModel2, {}); %we now need to figure out what happened to the essential rxns in the merge step above. %Note that it would be ideal to run checkTasks on the result from mergeLinear, but that doesn't work. @@ -339,6 +339,7 @@ prepData.minModel = scaledMinModel; prepData.refModelWithBM = closeModel(cModel); %do this here so we don't have to do this for each sample prepData.groupIds = groupIds; +prepData.reversedRxns = reversedRxns; %true for a reaction that runs reversed in its merged reaction prepData.essentialMetsForTasks = essentialMetsForTasks; prepData.toIgnoreExch = toIgnoreExch; From 221c815ff2b4353b78a7a74c1be1bc73fedacf56 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 29 Oct 2026 07:43:24 +0100 Subject: [PATCH 21/27] ftINIT: resolveTies orients a merged reaction as its smallest member and skips the flux phases in the last step --- CHANGELOG.md | 9 ++ src/INIT/ftINIT.m | 47 +++++--- src/INIT/ftINITInternalAlg.m | 210 ++++++++++++++++++++--------------- 3 files changed, 158 insertions(+), 108 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 525e2813..f6c0389c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -32,6 +32,15 @@ tie-break phases are skipped and the primary solution is kept, as in raven-toolbox (was: the phases ran on the unproven incumbent). The count and id-rank phases get half of `params.TimeLimit` each, rounded down and at least 1 s (was: the full limit). +* `ftINIT` with `resolveTies`: the direction phase minimises the backward flux of a merged + reaction in the orientation of its member with the smallest original id, and a turned-on + reaction with zero flux is forced in that orientation, so neither depends on which member + names the merged reaction. `prepINITModel` returns `reversedRxns` of `mergeLinear` in + `prepData`; with a `prepData` without it, every reaction keeps its own orientation. + `ftINITInternalAlg` takes the orientation as the new argument `rxnOrientation`. +* `ftINIT` with `resolveTies`: the last step skips the minimal-flux and direction phases, as + no step reads its fluxes (new `ftINITInternalAlg` argument `pinFluxes`). The reaction set + does not change. * `mergeLinear`: groups are numbered 1, 2, ... in the order of the first reaction of each group in `origRxnIds` (was: the order in which the reactions were merged). * `mergeLinear`: a reaction flipped after the first merge pass is the one marked in diff --git a/src/INIT/ftINIT.m b/src/INIT/ftINIT.m index f8449204..526354cd 100644 --- a/src/INIT/ftINIT.m +++ b/src/INIT/ftINIT.m @@ -71,17 +71,19 @@ % if true, pin each step's degenerate MILP to a canonical answer instead % of relying on the solver Seed alone: hold the score objective at its % optimum, then minimise the count of removable ("on") reactions, then, -% among the sparsest, minimise their summed reaction-id rank; in each -% main-extraction step (ftINITInternalAlg) the fluxes, which set the -% direction and magnitude a reaction is forced at in the next step, then -% come from the least-total-flux solution for that on/off state and, -% among those, from the one with the least id-rank-weighted backward -% flux (a merged reaction ranked by the smallest original id in its -% group). Applies to every main-extraction step and to the final -% gap-filling step (fitTasks' resolveTies) -- see raven-gecko-parity#104 -% and raven-toolbox's resolve_ties. A step whose primary solve stops at -% the time limit keeps its primary solution without these phases; the -% count and id-rank phases get half of the step's TimeLimit each, and a warning +% among the sparsest, minimise their summed reaction-id rank. This +% applies to every main-extraction step and to the final gap-filling +% step (fitTasks' resolveTies) -- see raven-gecko-parity#104 and +% raven-toolbox's resolve_ties. In each main-extraction step but the +% last (ftINITInternalAlg) the fluxes, which set the direction and +% magnitude a reaction is forced at in the next step, then come from the +% least-total-flux solution for that on/off state and, among those, from +% the one with the least id-rank-weighted backward flux (a merged +% reaction ranked by the smallest original id in its group and oriented +% as that member). The last step skips these flux phases, as no step +% reads its fluxes. A step whose primary solve stops at the time limit +% keeps its primary solution without any of these phases; the count and +% id-rank phases get half of the step's TimeLimit each, and a warning % names each phase that does not prove optimal or finds no solution % (default false). % proveAbsGap : double @@ -201,14 +203,23 @@ fluxes = ones(length(prepData.minModel.rxns), 1).*0.1; %resolveTies ranks a merged reaction by its group's smallest original id, -%which does not depend on which member names the merged reaction. +%which does not depend on which member names the merged reaction. The +%reference orientation of a merged reaction is that of the same member: +%-1 when that member runs reversed in the merged reaction. A prepData +%without reversedRxns gives every reaction orientation 1. rankKeys = prepData.minModel.rxns; +rxnOrientation = ones(numel(prepData.minModel.rxns),1); +hasReversed = isfield(prepData,'reversedRxns'); [inRef, refInd] = ismember(rankKeys, prepData.refModel.rxns); for i = find(inRef)' g = prepData.groupIds(refInd(i)); if g ~= 0 - members = sort(prepData.refModel.rxns(prepData.groupIds == g)); + memberInd = find(prepData.groupIds == g); + [members, ord] = sort(prepData.refModel.rxns(memberInd)); rankKeys{i} = members{1}; + if hasReversed && prepData.reversedRxns(memberInd(ord(1))) + rxnOrientation(i) = -1; + end end end @@ -220,6 +231,9 @@ error('RAVEN:badInput', '%s', 'RxnsToIgnoreMask may not cover rxns not covered in previous steps, but the other way around is fine.'); end rxnsToIgnoreLastStep = stp.RxnsToIgnoreMask; + %The fluxes of a step set the forced directions of the next one, so the + %last step skips the flux phases of resolveTies. + pinFluxes = initStep < numel(INITSteps); mm = prepData.minModel; @@ -249,8 +263,9 @@ %in the direction that they were previously carrying flux %first reverse the reactions that need to be reversed + %A zero flux takes the reference orientation of the reaction. rev = mm.rev == 1; - toRev = rxnsTurnedOn & rev & fluxes < 0; + toRev = rxnsTurnedOn & rev & (fluxes < 0 | (fluxes == 0 & rxnOrientation < 0)); mm = reverseRxns(mm, mm.rxns(toRev)); %Then make them irreversible @@ -277,7 +292,7 @@ params.MIPGap = 0; params.MIPGapAbs = proveAbsGap; try - [deletedRxnsInINIT1, metProduction,fullMipRes,rxnsTurnedOn1,fluxes1] = ftINITInternalAlg(mm,rxnScores,metData,essentialRxns,5,stp.AllowMetSecr,stp.PosRevOff,params, [], fluxes, verbose, 0.1, resolveTies, rankKeys); + [deletedRxnsInINIT1, metProduction,fullMipRes,rxnsTurnedOn1,fluxes1] = ftINITInternalAlg(mm,rxnScores,metData,essentialRxns,5,stp.AllowMetSecr,stp.PosRevOff,params, [], fluxes, verbose, 0.1, resolveTies, rankKeys, rxnOrientation, pinFluxes); fluxes1(toRev) = -fluxes1(toRev); success = true; if fullMipRes.stat ~= 1 @@ -331,7 +346,7 @@ if ~isempty(fullMipRes) startVals = fullMipRes.full; end - [deletedRxnsInINIT1, metProduction,fullMipRes,rxnsTurnedOn1,fluxes1] = ftINITInternalAlg(mm,rxnScores,metData,essentialRxns,5,stp.AllowMetSecr,stp.PosRevOff,params, startVals, fluxes, verbose, 0.1, resolveTies, rankKeys); + [deletedRxnsInINIT1, metProduction,fullMipRes,rxnsTurnedOn1,fluxes1] = ftINITInternalAlg(mm,rxnScores,metData,essentialRxns,5,stp.AllowMetSecr,stp.PosRevOff,params, startVals, fluxes, verbose, 0.1, resolveTies, rankKeys, rxnOrientation, pinFluxes); %This is a bit tricky - since we reversed some reactions, those fluxes also need to be reversed fluxes1(toRev) = -fluxes1(toRev); diff --git a/src/INIT/ftINITInternalAlg.m b/src/INIT/ftINITInternalAlg.m index 6caccf0e..99570503 100644 --- a/src/INIT/ftINITInternalAlg.m +++ b/src/INIT/ftINITInternalAlg.m @@ -1,4 +1,4 @@ -function [deletedRxns,metProduction,res,turnedOnRxns,fluxes]=ftINITInternalAlg(model,rxnScores,metData,essentialRxns,prodWeight,allowExcretion,remPosRev,params,startVals,fluxes,verbose,forceOn,resolveTies,rankKeys) +function [deletedRxns,metProduction,res,turnedOnRxns,fluxes]=ftINITInternalAlg(model,rxnScores,metData,essentialRxns,prodWeight,allowExcretion,remPosRev,params,startVals,fluxes,verbose,forceOn,resolveTies,rankKeys,rxnOrientation,pinFluxes) % ftINITInternalAlg % This function runs the MILP for a step in ftINIT. % @@ -63,6 +63,18 @@ % (optional, default model.rxns). ftINIT passes the smallest % original reaction id of each merge group, so the order % does not depend on which member names the merged reaction +% rxnOrientation vector with, for each reaction in model.rxns, 1 or -1: +% the sign that turns its flux into the direction of its +% reference orientation. Phase 2d of resolveTies minimises +% the backward flux in that orientation (optional, default +% all 1). ftINIT passes -1 for a merged reaction in which +% the member with the smallest original id runs reversed, +% so the preferred direction does not depend on which member +% names the merged reaction +% pinFluxes if false, resolveTies runs phases 2a and 2b only and the +% fluxes are those of phase 2b (optional, default true). +% ftINIT passes false in its last step, whose fluxes no +% later step reads % % deletedRxns reactions which were deleted by the algorithm (only % rxns included in the problem) @@ -92,6 +104,13 @@ if nargin < 13 || isempty(resolveTies) resolveTies = false; end +if nargin < 15 || isempty(rxnOrientation) + rxnOrientation = ones(numel(model.rxns),1); +end +rxnOrientation = rxnOrientation(:); +if nargin < 16 || isempty(pinFluxes) + pinFluxes = true; +end if isempty(essentialRxns) essentialRxns={}; end @@ -549,97 +568,104 @@ %Adopt the tie-resolved solution for the on/off read-out below. res = res2b; - %Phase 2c: the least total absolute flux for this on/off state. - %The next step forces each turned-on reaction in the direction - %of its flux, at min(0.99*|flux|,0.1), so the flux read-out - %comes from a minimal-flux solution instead of an arbitrary - %vertex. A binary on/off variable is fixed at its value; a - %continuous one that is on keeps at least its value, one that - %is off at most its value. Directions of reversible scored - %reactions stay free, so this is still a MILP. - %It stops at a relative or absolute gap of 1e-6, as raven-toolbox's - %phase 2c. - nCols = numel(prob.c); - minProb = tieProb2b; - if isfield(minProb,'start') - minProb = rmfield(minProb,'start'); - end - x = res2b.full; - negVal = round(x(negCols)); - minProb.lb(negCols) = negVal; - minProb.ub(negCols) = negVal; - posCols = [onoffPosIrrev, onoffPosRev]; - posVal = min(max(x(posCols),0),1); - posOn = posVal >= 0.5; - minProb.lb(posCols(posOn)) = posVal(posOn); - minProb.ub(posCols(~posOn)) = posVal(~posOn); - %|v| is v for a forward-only flux, -v for a backward-only one, - %and an added t >= v, t >= -v otherwise. - fwdOnly = minProb.lb(1:nRxns) >= 0; - bwdOnly = ~fwdOnly & minProb.ub(1:nRxns) <= 0; - tRxns = find(~fwdOnly & ~bwdOnly); - nT = numel(tRxns); - selV = sparse(1:nT, tRxns, 1, nT, nCols); - minProb.A = [minProb.A sparse(size(minProb.A,1),nT); -selV speye(nT); selV speye(nT)]; - minProb.a = minProb.A; - minProb.b = [minProb.b; zeros(2*nT,1)]; - minProb.csense = [minProb.csense repmat('G',1,2*nT)]; - minProb.lb = [minProb.lb; zeros(nT,1)]; - minProb.ub = [minProb.ub; inf(nT,1)]; - minProb.vartype = [minProb.vartype repmat('C',1,nT)]; - minProb.c = [zeros(nCols,1); ones(nT,1)]; - minProb.c(fwdOnly) = 1; - minProb.c(bwdOnly) = -1; - minParams = params; - minParams.MIPGap = 1e-6; - minParams.MIPGapAbs = 1e-6; - res2c = optimizeProb(minProb,minParams,verbose); - [isFeasible2c, isOptimal2c] = checkSolution(res2c); - if ~isOptimal2c - unproven{end+1} = 'phase2c-minflux'; - end - if isFeasible2c - %Only the fluxes change; the status stays that of phase 2b. - res.full = res2c.full(1:nCols); - - %Phase 2d: the least total flux can leave a direction open - %when two directions carry the same total flux, for example - %a pair of reversible reactions closing a loop either way. - %Keep the on/off state of phase 2c, hold the total absolute - %flux at most at its phase-2c value plus - %max(1e-7*total, 1e-6), and minimise the backward flux - %max(0,-v) of the reactions whose flux can take either sign, - %weighted by id rank from the end: sorted by rankKeys, the - %first of n reactions gets weight n and the last 1, so a - %reaction with a lower id prefers its forward direction. - if nT > 0 - total = minProb.c' * res2c.full; - dirProb = minProb; - nCols2c = numel(minProb.c); - [~, ord] = sort(rankKeys(tRxns)); - w = zeros(nT,1); - w(ord) = (nT:-1:1)'; - dirProb.A = [minProb.A sparse(size(minProb.A,1),nT); ... - minProb.c' sparse(1,nT); ... - selV sparse(nT,nT) speye(nT)]; - dirProb.a = dirProb.A; - dirProb.b = [minProb.b; total + max(abs(total)*1e-7, 1e-6); zeros(nT,1)]; - dirProb.csense = [minProb.csense 'L' repmat('G',1,nT)]; - dirProb.lb = [minProb.lb; zeros(nT,1)]; - dirProb.ub = [minProb.ub; inf(nT,1)]; - dirProb.vartype = [minProb.vartype repmat('C',1,nT)]; - dirProb.c = [zeros(nCols2c,1); w]; - %The phase-2c solution is the MIP start: the tight cap on - %the total flux leaves the solver's heuristics little room - %to find a first solution. - dirProb.start = [res2c.full; max(0, -res2c.full(tRxns))]; - res2d = optimizeProb(dirProb,minParams,verbose); - [isFeasible2d, isOptimal2d] = checkSolution(res2d); - if ~isOptimal2d - unproven{end+1} = 'phase2d-direction'; - end - if isFeasible2d - res.full = res2d.full(1:nCols); + %Phases 2c and 2d set the fluxes only, which the next ftINIT step + %reads; pinFluxes is false when no step follows. + if pinFluxes + %Phase 2c: the least total absolute flux for this on/off state. + %The next step forces each turned-on reaction in the direction + %of its flux, at min(0.99*|flux|,0.1), so the flux read-out + %comes from a minimal-flux solution instead of an arbitrary + %vertex. A binary on/off variable is fixed at its value; a + %continuous one that is on keeps at least its value, one that + %is off at most its value. Directions of reversible scored + %reactions stay free, so this is still a MILP. + %It stops at a relative or absolute gap of 1e-6, as raven-toolbox's + %phase 2c. + nCols = numel(prob.c); + minProb = tieProb2b; + if isfield(minProb,'start') + minProb = rmfield(minProb,'start'); + end + x = res2b.full; + negVal = round(x(negCols)); + minProb.lb(negCols) = negVal; + minProb.ub(negCols) = negVal; + posCols = [onoffPosIrrev, onoffPosRev]; + posVal = min(max(x(posCols),0),1); + posOn = posVal >= 0.5; + minProb.lb(posCols(posOn)) = posVal(posOn); + minProb.ub(posCols(~posOn)) = posVal(~posOn); + %|v| is v for a forward-only flux, -v for a backward-only one, + %and an added t >= v, t >= -v otherwise. + fwdOnly = minProb.lb(1:nRxns) >= 0; + bwdOnly = ~fwdOnly & minProb.ub(1:nRxns) <= 0; + tRxns = find(~fwdOnly & ~bwdOnly); + nT = numel(tRxns); + selV = sparse(1:nT, tRxns, 1, nT, nCols); + minProb.A = [minProb.A sparse(size(minProb.A,1),nT); -selV speye(nT); selV speye(nT)]; + minProb.a = minProb.A; + minProb.b = [minProb.b; zeros(2*nT,1)]; + minProb.csense = [minProb.csense repmat('G',1,2*nT)]; + minProb.lb = [minProb.lb; zeros(nT,1)]; + minProb.ub = [minProb.ub; inf(nT,1)]; + minProb.vartype = [minProb.vartype repmat('C',1,nT)]; + minProb.c = [zeros(nCols,1); ones(nT,1)]; + minProb.c(fwdOnly) = 1; + minProb.c(bwdOnly) = -1; + minParams = params; + minParams.MIPGap = 1e-6; + minParams.MIPGapAbs = 1e-6; + res2c = optimizeProb(minProb,minParams,verbose); + [isFeasible2c, isOptimal2c] = checkSolution(res2c); + if ~isOptimal2c + unproven{end+1} = 'phase2c-minflux'; + end + if isFeasible2c + %Only the fluxes change; the status stays that of phase 2b. + res.full = res2c.full(1:nCols); + + %Phase 2d: the least total flux can leave a direction open + %when two directions carry the same total flux, for example + %a pair of reversible reactions closing a loop either way. + %Keep the on/off state of phase 2c, hold the total absolute + %flux at most at its phase-2c value plus + %max(1e-7*total, 1e-6), and minimise the backward flux + %max(0,-o*v) of the reactions whose flux can take either + %sign, with o their rxnOrientation, weighted by id rank + %from the end: sorted by rankKeys, the first of n + %reactions gets weight n and the last 1, so a reaction + %with a lower id prefers its forward direction. + if nT > 0 + total = minProb.c' * res2c.full; + dirProb = minProb; + nCols2c = numel(minProb.c); + [~, ord] = sort(rankKeys(tRxns)); + w = zeros(nT,1); + w(ord) = (nT:-1:1)'; + o = rxnOrientation(tRxns); + selO = sparse(1:nT, tRxns, o, nT, nCols2c); + dirProb.A = [minProb.A sparse(size(minProb.A,1),nT); ... + minProb.c' sparse(1,nT); ... + selO speye(nT)]; + dirProb.a = dirProb.A; + dirProb.b = [minProb.b; total + max(abs(total)*1e-7, 1e-6); zeros(nT,1)]; + dirProb.csense = [minProb.csense 'L' repmat('G',1,nT)]; + dirProb.lb = [minProb.lb; zeros(nT,1)]; + dirProb.ub = [minProb.ub; inf(nT,1)]; + dirProb.vartype = [minProb.vartype repmat('C',1,nT)]; + dirProb.c = [zeros(nCols2c,1); w]; + %The phase-2c solution is the MIP start: the tight cap on + %the total flux leaves the solver's heuristics little room + %to find a first solution. + dirProb.start = [res2c.full; max(0, -o.*res2c.full(tRxns))]; + res2d = optimizeProb(dirProb,minParams,verbose); + [isFeasible2d, isOptimal2d] = checkSolution(res2d); + if ~isOptimal2d + unproven{end+1} = 'phase2d-direction'; + end + if isFeasible2d + res.full = res2d.full(1:nCols); + end end end end From ccf636104a132040f4df14f862806e9b285873c4 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 29 Oct 2026 07:43:24 +0100 Subject: [PATCH 22/27] tINIT: resolveTies follows rxnOrientation, pinFluxes skips the flux phases, prepData carries reversedRxns --- tests/function_tests/tINIT.m | 179 +++++++++++++++++++++++++++-------- 1 file changed, 139 insertions(+), 40 deletions(-) diff --git a/tests/function_tests/tINIT.m b/tests/function_tests/tINIT.m index 4c4e9c32..cf8b5496 100644 --- a/tests/function_tests/tINIT.m +++ b/tests/function_tests/tINIT.m @@ -520,18 +520,7 @@ function ftINITInternalAlgResolveTiesFixesEqualFluxDirections(testCase) % flux: R2 when ranked by id, R1 when rankKeys rank it last. N is % a negative-score dead end. testCase.assumeMILPSolver(); - m.rxns = {'IN_x';'IN_z';'R1';'R2';'E';'OUT_p';'N'}; - m.mets = {'x';'y';'z';'p';'d'}; - % INx INz R1 R2 E OUTp N - m.S = sparse([1 0 1 0 0 0 -1; % x - 0 0 -1 -1 -1 0 0; % y - 0 1 0 1 0 0 0; % z - 0 0 0 0 1 -1 0; % p - 0 0 0 0 0 0 1]); % d - m.lb = [0;0;-1000;-1000;0;0;0]; - m.ub = 1000*ones(7,1); - m.rev = double(m.lb < 0); - scores = [0;0;0;0;0;0;-1]; + [m, scores] = getTstModelTwoRoutes(); fluxIn = ones(7,1); [~,~,~,~,fluxes] = ftINITInternalAlg(m,scores,[],{'E'}, ... [],false,false,struct(),[],fluxIn,false,0.1,true); @@ -545,35 +534,12 @@ function ftINITInternalAlgResolveTiesFixesEqualFluxDirections(testCase) function ftINITInternalAlgResolveTiesSkippedAfterTimeLimit(testCase) % When the primary solve stops at the time limit, resolveTies % keeps the primary solution and runs no tie-break phase. A - % wrapper of optimizeProb earlier on the path counts the solves - % and reports each as stopped at the time limit. + % wrapper of optimizeProb counts the solves and reports each as + % stopped at the time limit. testCase.assumeMILPSolver(); - global TINIT_OPTIMIZEPROB TINIT_OPTIMIZEPROB_CALLS - TINIT_OPTIMIZEPROB = @optimizeProb; - TINIT_OPTIMIZEPROB_CALLS = 0; - d = tempname; mkdir(d); - fid = fopen(fullfile(d,'optimizeProb.m'),'w'); - fprintf(fid,'%s\n', ... - 'function res = optimizeProb(varargin)', ... - 'global TINIT_OPTIMIZEPROB TINIT_OPTIMIZEPROB_CALLS', ... - 'res = TINIT_OPTIMIZEPROB(varargin{:});', ... - 'TINIT_OPTIMIZEPROB_CALLS = TINIT_OPTIMIZEPROB_CALLS + 1;', ... - 'res.hitTimeLimit = true;', ... - 'end'); - fclose(fid); - testCase.addTeardown(@() rmdir(d,'s')); - addpath(d,'-begin'); - testCase.addTeardown(@() rmpath(d)); - testCase.addTeardown(@() clear('global','TINIT_OPTIMIZEPROB','TINIT_OPTIMIZEPROB_CALLS')); - % R1 and R2 can each feed the essential E: a tie for resolveTies. - m.rxns = {'R1';'R2';'E';'EX_a';'EX_p'}; - m.mets = {'a';'m';'p'}; - m.S = sparse([-1 -1 0 -1 0; - 1 1 -1 0 0; - 0 0 1 0 -1]); - m.lb = [0;0;0;-1000;-1000]; - m.ub = 1000*ones(5,1); - m.rev = double(m.lb < 0); + installOptimizeProbCounter(testCase, true); + global TINIT_OPTIMIZEPROB_CALLS + m = getTstModelTwoSources(); [~,~,res,turnedOn] = ftINITInternalAlg(m,[-1;-1;0;0;0],[],{'E'}, ... [],false,false,struct(),[],0.1*ones(5,1),false,0.1,true); testCase.verifyEqual(TINIT_OPTIMIZEPROB_CALLS, 1); @@ -581,6 +547,48 @@ function ftINITInternalAlgResolveTiesSkippedAfterTimeLimit(testCase) testCase.verifyEqual(nnz(turnedOn(1:2)), 1); end + function ftINITInternalAlgResolveTiesWithoutPinFluxesSkipsFluxPhases(testCase) + % With pinFluxes false, resolveTies runs the primary solve and + % phases 2a and 2b only, and still keeps the lower-id R1. + testCase.assumeMILPSolver(); + installOptimizeProbCounter(testCase, false); + global TINIT_OPTIMIZEPROB_CALLS + m = getTstModelTwoSources(); + [~,~,~,turnedOn] = ftINITInternalAlg(m,[-1;-1;0;0;0],[],{'E'}, ... + [],false,false,struct(),[],0.1*ones(5,1),false,0.1,true,[],[],false); + testCase.verifyEqual(TINIT_OPTIMIZEPROB_CALLS, 3); + testCase.verifyEqual(turnedOn(1:2), [true false]); + TINIT_OPTIMIZEPROB_CALLS = 0; + ftINITInternalAlg(m,[-1;-1;0;0;0],[],{'E'}, ... + [],false,false,struct(),[],0.1*ones(5,1),false,0.1,true); + testCase.verifyEqual(TINIT_OPTIMIZEPROB_CALLS, 5); + end + + function ftINITInternalAlgResolveTiesFollowsRxnOrientation(testCase) + % The two-route model of + % ftINITInternalAlgResolveTiesFixesEqualFluxDirections: R1 ranks + % first, so by default R2 carries the backward flux. With + % rxnOrientation -1 for R1, its backward flux in that orientation + % is a forward flux, so R1 runs backward instead. + testCase.assumeMILPSolver(); + [m, scores] = getTstModelTwoRoutes(); + orient = [1;1;-1;1;1;1;1]; + [~,~,~,~,fluxes] = ftINITInternalAlg(m,scores,[],{'E'}, ... + [],false,false,struct(),[],ones(7,1),false,0.1,true,[],orient); + testCase.verifyEqual(fluxes(3:4), [-0.1;0], 'AbsTol', 1e-6); + end + + function prepINITModelReturnsReversedRxnsOfTheMerge(testCase) + % Rz (h[c] => y[c]) and Ra (g[c] => y[c]) both produce y[c], so + % the merge runs Ra reversed in the merged reaction Rz; Rp and + % Rq form a second route between h[c] and g[c]. + testModel = getTstModelFlippedMember(); + evalc('prepData = prepINITModel(testModel, [], {}, false, {}, ''s'');'); + testCase.verifyClass(prepData.reversedRxns, 'logical'); + testCase.verifyEqual(prepData.refModel.rxns(prepData.reversedRxns), {'Ra'}); + testCase.verifyTrue(ismember('Rz', prepData.minModel.rxns)); + end + function ftINITWithTaskRuns(testCase) % A task requiring e[s] from a[s] forces R2 and R7 to be essential. testCase.assumeMILPSolver(); @@ -1013,6 +1021,97 @@ function scoreModelRejectsUnknownDataPrecedence(testCase) testModel.b = repmat(0,8,1); end +function installOptimizeProbCounter(testCase, fakeTimeLimit) +%Put a wrapper of optimizeProb first on the path for the rest of the test. It +%counts the solves in the global TINIT_OPTIMIZEPROB_CALLS and, if +%fakeTimeLimit, reports each solve as stopped at the time limit. +global TINIT_OPTIMIZEPROB TINIT_OPTIMIZEPROB_CALLS TINIT_OPTIMIZEPROB_FAKETL +TINIT_OPTIMIZEPROB = @optimizeProb; +TINIT_OPTIMIZEPROB_CALLS = 0; +TINIT_OPTIMIZEPROB_FAKETL = fakeTimeLimit; +d = tempname; mkdir(d); +fid = fopen(fullfile(d,'optimizeProb.m'),'w'); +fprintf(fid,'%s\n', ... + 'function res = optimizeProb(varargin)', ... + 'global TINIT_OPTIMIZEPROB TINIT_OPTIMIZEPROB_CALLS TINIT_OPTIMIZEPROB_FAKETL', ... + 'res = TINIT_OPTIMIZEPROB(varargin{:});', ... + 'TINIT_OPTIMIZEPROB_CALLS = TINIT_OPTIMIZEPROB_CALLS + 1;', ... + 'if TINIT_OPTIMIZEPROB_FAKETL, res.hitTimeLimit = true; end', ... + 'end'); +fclose(fid); +testCase.addTeardown(@() rmdir(d,'s')); +addpath(d,'-begin'); +testCase.addTeardown(@() rmpath(d)); +testCase.addTeardown(@() clear('global','TINIT_OPTIMIZEPROB','TINIT_OPTIMIZEPROB_CALLS','TINIT_OPTIMIZEPROB_FAKETL')); +end + +function m = getTstModelTwoSources() +%R1 and R2 (a => m) can each feed the essential E (m => p): a tie for +%resolveTies. +m.rxns = {'R1';'R2';'E';'EX_a';'EX_p'}; +m.mets = {'a';'m';'p'}; +m.S = sparse([-1 -1 0 -1 0; + 1 1 -1 0 0; + 0 0 1 0 -1]); +m.lb = [0;0;0;-1000;-1000]; +m.ub = 1000*ones(5,1); +m.rev = double(m.lb < 0); +end + +function [m, scores] = getTstModelTwoRoutes() +%The essential E consumes y, which the reversible R1 (y <=> x) or R2 +%(y <=> z) can supply when run backward from the imports IN_x or IN_z. N is a +%negative-score dead end. +m.rxns = {'IN_x';'IN_z';'R1';'R2';'E';'OUT_p';'N'}; +m.mets = {'x';'y';'z';'p';'d'}; +% INx INz R1 R2 E OUTp N +m.S = sparse([1 0 1 0 0 0 -1; % x + 0 0 -1 -1 -1 0 0; % y + 0 1 0 1 0 0 0; % z + 0 0 0 0 1 -1 0; % p + 0 0 0 0 0 0 1]); % d +m.lb = [0;0;-1000;-1000;0;0;0]; +m.ub = 1000*ones(7,1); +m.rev = double(m.lb < 0); +scores = [0;0;0;0;0;0;-1]; +end + +function testModel = getTstModelFlippedMember() +%h[s] and g[s] are exchanged; h[c] reaches g[c] through y[c] (Rz, Ra) or +%through y2[c] (Rp, Rq). Rz and Ra are both written as producers of y[c]. +testModel = struct(); +testModel.id = 'flippedMember'; +testModel.rxns = {}; +testModel.S = []; +testModel.rev = []; +testModel.mets = {'hs';'hc';'gc';'y1';'y2';'gs'}; +testModel.metNames = {'h';'h';'g';'y';'y2';'g'}; +testModel.comps = {'s';'c'}; +testModel.compNames = testModel.comps; +testModel.metComps = [1;2;2;2;2;1]; +testModel.genes = {'Gz';'Ga';'Gp';'Gq'}; +testModel.grRules = {}; +testModel.rxnGeneMat = []; + +rxnsToAdd = struct(); +rxnsToAdd.rxns = {'EX_h';'Th';'Rz';'Ra';'Rp';'Rq';'Tg';'EX_g'}; +rxnsToAdd.equations = {'<=> h[s]';... + 'h[s] <=> h[c]';... + 'h[c] <=> y[c]';... + 'g[c] <=> y[c]';... + 'h[c] <=> y2[c]';... + 'y2[c] <=> g[c]';... + 'g[c] <=> g[s]';... + 'g[s] <=>'}; +rxnsToAdd.grRules = {'';'';'Gz';'Ga';'Gp';'Gq';'';''}; +evalc('testModel = addRxns(testModel,rxnsToAdd, 3);'); +testModel.c = [0;0;0;0;0;0;0;1]; +testModel.ub = repmat(1000,8,1); +testModel.lb = repmat(-1000,8,1); +testModel.rxnNames = testModel.rxns; +testModel.b = zeros(6,1); +end + function testModel = getTstModelTwoPaths() %a[c] reaches m[c] by Rm, Rc or by Rn, Rb; E, R7 and R8 export e[s]. The %paths are reversible: simplifyModel in prepINITModel fails on a model From 8307a1526376d629bbaef83ff570120e93929d11 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 29 Oct 2026 07:50:02 +0100 Subject: [PATCH 23/27] simplifyModel: constrainReversible changes nothing on a model without a reversible reaction --- CHANGELOG.md | 5 ++++ src/manipulation/simplifyModel.m | 41 +++++++++++++++++--------------- 2 files changed, 27 insertions(+), 19 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index f6c0389c..1a0f8ebe 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,11 @@ * `mergeCompartments`: a merged metabolite is the copy already in compartment `s`, if any, otherwise the copy whose id sorts first, with every field of that copy (was: the first copy in model order), so `ACCOA_c` and `ACCOA_m` become `ACCOA_c` whatever their order. +* `simplifyModel` with `constrainReversible`: on a model without a reversible reaction this step + changes nothing (was: an error "logical indices ... outside of the array bounds", as + `getAllowedBounds` took the empty reaction list for all reactions). `prepINITModel` therefore + runs on a model whose reactions are all irreversible after the first simplification, with the + result of raven-toolbox's `prep_init_model`. * `removeGenes`: the help states that, when at least one of the requested genes is in the model, genes no remaining reaction uses are removed as well; when none is, no gene, reaction or bound changes. diff --git a/src/manipulation/simplifyModel.m b/src/manipulation/simplifyModel.m index e5dd214f..d579230d 100755 --- a/src/manipulation/simplifyModel.m +++ b/src/manipulation/simplifyModel.m @@ -185,30 +185,33 @@ if constrainReversible==true revs=find(reducedModel.rev); - [I,J]=getAllowedBounds(reducedModel,revs); + %getAllowedBounds treats an empty reaction list as all reactions + if ~isempty(revs) + [I,J]=getAllowedBounds(reducedModel,revs); - I=abs(I); - J=abs(J); + I=abs(I); + J=abs(J); - %Get the "small" values - K=I<10^-10; - L=J<10^-10; + %Get the "small" values + K=I<10^-10; + L=J<10^-10; - %Keep the values where only one direction is zero (small) - I(K==L)=[]; - J(K==L)=[]; - revs(K==L)=[]; + %Keep the values where only one direction is zero (small) + I(K==L)=[]; + J(K==L)=[]; + revs(K==L)=[]; - %Change the reversibility of the remaining reactions - reducedModel.rev(revs(J>10^-10))=0; %Ignore reverse direction - reducedModel.lb(revs(J>10^-10))=0; + %Change the reversibility of the remaining reactions + reducedModel.rev(revs(J>10^-10))=0; %Ignore reverse direction + reducedModel.lb(revs(J>10^-10))=0; - toSwitch=revs(I>10^-10); - reducedModel.rev(toSwitch)=0; %Change directionality - reducedModel.ub(toSwitch)=reducedModel.lb(toSwitch)*-1; - reducedModel.lb(toSwitch)=0; - reducedModel.S(:,toSwitch)=reducedModel.S(:,toSwitch).*-1; - reducedModel.c(toSwitch)=reducedModel.c(toSwitch)*-1; + toSwitch=revs(I>10^-10); + reducedModel.rev(toSwitch)=0; %Change directionality + reducedModel.ub(toSwitch)=reducedModel.lb(toSwitch)*-1; + reducedModel.lb(toSwitch)=0; + reducedModel.S(:,toSwitch)=reducedModel.S(:,toSwitch).*-1; + reducedModel.c(toSwitch)=reducedModel.c(toSwitch)*-1; + end end if groupLinear==true if suppressWarnings==false From 6926a66dbbed8984fc902dcccc76b5c69d50ab67 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 29 Oct 2026 07:50:02 +0100 Subject: [PATCH 24/27] tINIT: prepINITModel runs on a model without a reversible reaction --- tests/function_tests/tINIT.m | 49 ++++++++++++++++++++++++++++++++++-- 1 file changed, 47 insertions(+), 2 deletions(-) diff --git a/tests/function_tests/tINIT.m b/tests/function_tests/tINIT.m index cf8b5496..a89beef1 100644 --- a/tests/function_tests/tINIT.m +++ b/tests/function_tests/tINIT.m @@ -578,6 +578,17 @@ function ftINITInternalAlgResolveTiesFollowsRxnOrientation(testCase) testCase.verifyEqual(fluxes(3:4), [-0.1;0], 'AbsTol', 1e-6); end + function prepINITModelHandlesModelWithoutReversibleRxns(testCase) + % Every reaction is irreversible, so the second simplification + % has no reversible reaction to constrain. The merge gives the + % same reactions and groups as raven-toolbox's prep_init_model. + testModel = getTstModelIrreversible(); + evalc('prepData = prepINITModel(testModel, [], {}, false, {}, ''s'');'); + testCase.verifyEqual(prepData.minModel.rxns, {'EX_a';'R1';'R2';'T2'}); + testCase.verifyEqual(prepData.groupIds, [1;1;0;2;2;3;3]); + testCase.verifyEqual(prepData.refModel.rxns, testModel.rxns); + end + function prepINITModelReturnsReversedRxnsOfTheMerge(testCase) % Rz (h[c] => y[c]) and Ra (g[c] => y[c]) both produce y[c], so % the merge runs Ra reversed in the merged reaction Rz; Rp and @@ -1076,6 +1087,41 @@ function installOptimizeProbCounter(testCase, fakeTimeLimit) scores = [0;0;0;0;0;0;-1]; end +function testModel = getTstModelIrreversible() +%a[s] is imported and reaches b[c] by R1 or by R2, R3; b[c] is exported. +%Every reaction is irreversible. +testModel = struct(); +testModel.id = 'irreversible'; +testModel.rxns = {}; +testModel.S = []; +testModel.rev = []; +testModel.mets = {'as';'ac';'bc';'cc';'bs'}; +testModel.metNames = {'a';'a';'b';'c';'b'}; +testModel.comps = {'s';'c'}; +testModel.compNames = testModel.comps; +testModel.metComps = [1;2;2;2;1]; +testModel.genes = {'G1';'G2';'G3'}; +testModel.grRules = {}; +testModel.rxnGeneMat = []; + +rxnsToAdd = struct(); +rxnsToAdd.rxns = {'EX_a';'T1';'R1';'R2';'R3';'T2';'EX_b'}; +rxnsToAdd.equations = {'=> a[s]';... + 'a[s] => a[c]';... + 'a[c] => b[c]';... + 'a[c] => c[c]';... + 'c[c] => b[c]';... + 'b[c] => b[s]';... + 'b[s] =>'}; +rxnsToAdd.grRules = {'';'';'G1';'G2';'G3';'';''}; +evalc('testModel = addRxns(testModel,rxnsToAdd, 3);'); +testModel.c = [0;0;0;0;0;0;1]; +testModel.ub = repmat(1000,7,1); +testModel.lb = zeros(7,1); +testModel.rxnNames = testModel.rxns; +testModel.b = zeros(5,1); +end + function testModel = getTstModelFlippedMember() %h[s] and g[s] are exchanged; h[c] reaches g[c] through y[c] (Rz, Ra) or %through y2[c] (Rp, Rq). Rz and Ra are both written as producers of y[c]. @@ -1114,8 +1160,7 @@ function installOptimizeProbCounter(testCase, fakeTimeLimit) function testModel = getTstModelTwoPaths() %a[c] reaches m[c] by Rm, Rc or by Rn, Rb; E, R7 and R8 export e[s]. The -%paths are reversible: simplifyModel in prepINITModel fails on a model -%without a reversible reaction after the first simplification. +%paths are reversible. testModel = struct(); testModel.id = 'twoPaths'; testModel.rxns = {}; From a6dd6651e98f8235e6787e8c166989eab637caac Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 29 Oct 2026 10:22:32 +0100 Subject: [PATCH 25/27] ftINITInternalAlg: help and comments of pinFluxes name the phases it skips --- src/INIT/ftINITInternalAlg.m | 7 ++++--- tests/function_tests/tINIT.m | 2 +- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/src/INIT/ftINITInternalAlg.m b/src/INIT/ftINITInternalAlg.m index 99570503..1f2a9dd5 100644 --- a/src/INIT/ftINITInternalAlg.m +++ b/src/INIT/ftINITInternalAlg.m @@ -71,7 +71,7 @@ % the member with the smallest original id runs reversed, % so the preferred direction does not depend on which member % names the merged reaction -% pinFluxes if false, resolveTies runs phases 2a and 2b only and the +% pinFluxes if false, resolveTies skips phases 2c and 2d and the % fluxes are those of phase 2b (optional, default true). % ftINIT passes false in its last step, whose fluxes no % later step reads @@ -568,8 +568,9 @@ %Adopt the tie-resolved solution for the on/off read-out below. res = res2b; - %Phases 2c and 2d set the fluxes only, which the next ftINIT step - %reads; pinFluxes is false when no step follows. + %Phases 2c and 2d change the fluxes, not the on/off state; the next + %ftINIT step reads the fluxes, and pinFluxes is false when no step + %follows. if pinFluxes %Phase 2c: the least total absolute flux for this on/off state. %The next step forces each turned-on reaction in the direction diff --git a/tests/function_tests/tINIT.m b/tests/function_tests/tINIT.m index a89beef1..9d22b868 100644 --- a/tests/function_tests/tINIT.m +++ b/tests/function_tests/tINIT.m @@ -549,7 +549,7 @@ function ftINITInternalAlgResolveTiesSkippedAfterTimeLimit(testCase) function ftINITInternalAlgResolveTiesWithoutPinFluxesSkipsFluxPhases(testCase) % With pinFluxes false, resolveTies runs the primary solve and - % phases 2a and 2b only, and still keeps the lower-id R1. + % phases 2a and 2b, not 2c and 2d, and still keeps the lower-id R1. testCase.assumeMILPSolver(); installOptimizeProbCounter(testCase, false); global TINIT_OPTIMIZEPROB_CALLS From 1461f74148ec1a8158e011bcd512d90b22378bb3 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 29 Oct 2026 17:40:34 +0100 Subject: [PATCH 26/27] ftINIT: resolveTies keeps the orientation of each merged reaction as written Phase 2d of ftINITInternalAlg weights the backward flux by the rank key only, a turned-on reaction with zero flux is forced forward in the next step, and prepData has no reversedRxns field. ftINITInternalAlg takes pinFluxes as its 15th argument. tINIT checks that ftINIT passes pinFluxes false in the last step only. --- CHANGELOG.md | 9 +--- src/INIT/ftINIT.m | 30 ++++------- src/INIT/ftINITInternalAlg.m | 34 ++++-------- src/INIT/prepINITModel.m | 3 +- tests/function_tests/tINIT.m | 100 +++++++++++++---------------------- 5 files changed, 61 insertions(+), 115 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 1a0f8ebe..d047735b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -37,15 +37,10 @@ tie-break phases are skipped and the primary solution is kept, as in raven-toolbox (was: the phases ran on the unproven incumbent). The count and id-rank phases get half of `params.TimeLimit` each, rounded down and at least 1 s (was: the full limit). -* `ftINIT` with `resolveTies`: the direction phase minimises the backward flux of a merged - reaction in the orientation of its member with the smallest original id, and a turned-on - reaction with zero flux is forced in that orientation, so neither depends on which member - names the merged reaction. `prepINITModel` returns `reversedRxns` of `mergeLinear` in - `prepData`; with a `prepData` without it, every reaction keeps its own orientation. - `ftINITInternalAlg` takes the orientation as the new argument `rxnOrientation`. * `ftINIT` with `resolveTies`: the last step skips the minimal-flux and direction phases, as no step reads its fluxes (new `ftINITInternalAlg` argument `pinFluxes`). The reaction set - does not change. + does not change; the fluxes and `metProduction` of the last step are those of the id-rank + phase. * `mergeLinear`: groups are numbered 1, 2, ... in the order of the first reaction of each group in `origRxnIds` (was: the order in which the reactions were merged). * `mergeLinear`: a reaction flipped after the first merge pass is the one marked in diff --git a/src/INIT/ftINIT.m b/src/INIT/ftINIT.m index 526354cd..61de0623 100644 --- a/src/INIT/ftINIT.m +++ b/src/INIT/ftINIT.m @@ -79,11 +79,11 @@ % magnitude a reaction is forced at in the next step, then come from the % least-total-flux solution for that on/off state and, among those, from % the one with the least id-rank-weighted backward flux (a merged -% reaction ranked by the smallest original id in its group and oriented -% as that member). The last step skips these flux phases, as no step -% reads its fluxes. A step whose primary solve stops at the time limit -% keeps its primary solution without any of these phases; the count and -% id-rank phases get half of the step's TimeLimit each, and a warning +% reaction ranked by the smallest original id in its group). The last +% step skips these flux phases, as no step reads its fluxes. A step +% whose primary solve stops at the time limit keeps its primary +% solution without any of these phases; the count and id-rank phases +% get half of the step's TimeLimit each, and a warning % names each phase that does not prove optimal or finds no solution % (default false). % proveAbsGap : double @@ -203,23 +203,14 @@ fluxes = ones(length(prepData.minModel.rxns), 1).*0.1; %resolveTies ranks a merged reaction by its group's smallest original id, -%which does not depend on which member names the merged reaction. The -%reference orientation of a merged reaction is that of the same member: -%-1 when that member runs reversed in the merged reaction. A prepData -%without reversedRxns gives every reaction orientation 1. +%which does not depend on which member names the merged reaction. rankKeys = prepData.minModel.rxns; -rxnOrientation = ones(numel(prepData.minModel.rxns),1); -hasReversed = isfield(prepData,'reversedRxns'); [inRef, refInd] = ismember(rankKeys, prepData.refModel.rxns); for i = find(inRef)' g = prepData.groupIds(refInd(i)); if g ~= 0 - memberInd = find(prepData.groupIds == g); - [members, ord] = sort(prepData.refModel.rxns(memberInd)); + members = sort(prepData.refModel.rxns(prepData.groupIds == g)); rankKeys{i} = members{1}; - if hasReversed && prepData.reversedRxns(memberInd(ord(1))) - rxnOrientation(i) = -1; - end end end @@ -263,9 +254,8 @@ %in the direction that they were previously carrying flux %first reverse the reactions that need to be reversed - %A zero flux takes the reference orientation of the reaction. rev = mm.rev == 1; - toRev = rxnsTurnedOn & rev & (fluxes < 0 | (fluxes == 0 & rxnOrientation < 0)); + toRev = rxnsTurnedOn & rev & fluxes < 0; mm = reverseRxns(mm, mm.rxns(toRev)); %Then make them irreversible @@ -292,7 +282,7 @@ params.MIPGap = 0; params.MIPGapAbs = proveAbsGap; try - [deletedRxnsInINIT1, metProduction,fullMipRes,rxnsTurnedOn1,fluxes1] = ftINITInternalAlg(mm,rxnScores,metData,essentialRxns,5,stp.AllowMetSecr,stp.PosRevOff,params, [], fluxes, verbose, 0.1, resolveTies, rankKeys, rxnOrientation, pinFluxes); + [deletedRxnsInINIT1, metProduction,fullMipRes,rxnsTurnedOn1,fluxes1] = ftINITInternalAlg(mm,rxnScores,metData,essentialRxns,5,stp.AllowMetSecr,stp.PosRevOff,params, [], fluxes, verbose, 0.1, resolveTies, rankKeys, pinFluxes); fluxes1(toRev) = -fluxes1(toRev); success = true; if fullMipRes.stat ~= 1 @@ -346,7 +336,7 @@ if ~isempty(fullMipRes) startVals = fullMipRes.full; end - [deletedRxnsInINIT1, metProduction,fullMipRes,rxnsTurnedOn1,fluxes1] = ftINITInternalAlg(mm,rxnScores,metData,essentialRxns,5,stp.AllowMetSecr,stp.PosRevOff,params, startVals, fluxes, verbose, 0.1, resolveTies, rankKeys, rxnOrientation, pinFluxes); + [deletedRxnsInINIT1, metProduction,fullMipRes,rxnsTurnedOn1,fluxes1] = ftINITInternalAlg(mm,rxnScores,metData,essentialRxns,5,stp.AllowMetSecr,stp.PosRevOff,params, startVals, fluxes, verbose, 0.1, resolveTies, rankKeys, pinFluxes); %This is a bit tricky - since we reversed some reactions, those fluxes also need to be reversed fluxes1(toRev) = -fluxes1(toRev); diff --git a/src/INIT/ftINITInternalAlg.m b/src/INIT/ftINITInternalAlg.m index 1f2a9dd5..af29c682 100644 --- a/src/INIT/ftINITInternalAlg.m +++ b/src/INIT/ftINITInternalAlg.m @@ -1,4 +1,4 @@ -function [deletedRxns,metProduction,res,turnedOnRxns,fluxes]=ftINITInternalAlg(model,rxnScores,metData,essentialRxns,prodWeight,allowExcretion,remPosRev,params,startVals,fluxes,verbose,forceOn,resolveTies,rankKeys,rxnOrientation,pinFluxes) +function [deletedRxns,metProduction,res,turnedOnRxns,fluxes]=ftINITInternalAlg(model,rxnScores,metData,essentialRxns,prodWeight,allowExcretion,remPosRev,params,startVals,fluxes,verbose,forceOn,resolveTies,rankKeys,pinFluxes) % ftINITInternalAlg % This function runs the MILP for a step in ftINIT. % @@ -63,16 +63,9 @@ % (optional, default model.rxns). ftINIT passes the smallest % original reaction id of each merge group, so the order % does not depend on which member names the merged reaction -% rxnOrientation vector with, for each reaction in model.rxns, 1 or -1: -% the sign that turns its flux into the direction of its -% reference orientation. Phase 2d of resolveTies minimises -% the backward flux in that orientation (optional, default -% all 1). ftINIT passes -1 for a merged reaction in which -% the member with the smallest original id runs reversed, -% so the preferred direction does not depend on which member -% names the merged reaction % pinFluxes if false, resolveTies skips phases 2c and 2d and the -% fluxes are those of phase 2b (optional, default true). +% fluxes and metProduction are those of phase 2b +% (optional, default true). % ftINIT passes false in its last step, whose fluxes no % later step reads % @@ -104,11 +97,7 @@ if nargin < 13 || isempty(resolveTies) resolveTies = false; end -if nargin < 15 || isempty(rxnOrientation) - rxnOrientation = ones(numel(model.rxns),1); -end -rxnOrientation = rxnOrientation(:); -if nargin < 16 || isempty(pinFluxes) +if nargin < 15 || isempty(pinFluxes) pinFluxes = true; end if isempty(essentialRxns) @@ -631,11 +620,10 @@ %Keep the on/off state of phase 2c, hold the total absolute %flux at most at its phase-2c value plus %max(1e-7*total, 1e-6), and minimise the backward flux - %max(0,-o*v) of the reactions whose flux can take either - %sign, with o their rxnOrientation, weighted by id rank - %from the end: sorted by rankKeys, the first of n - %reactions gets weight n and the last 1, so a reaction - %with a lower id prefers its forward direction. + %max(0,-v) of the reactions whose flux can take either sign, + %weighted by id rank from the end: sorted by rankKeys, the + %first of n reactions gets weight n and the last 1, so a + %reaction with a lower id prefers its forward direction. if nT > 0 total = minProb.c' * res2c.full; dirProb = minProb; @@ -643,11 +631,9 @@ [~, ord] = sort(rankKeys(tRxns)); w = zeros(nT,1); w(ord) = (nT:-1:1)'; - o = rxnOrientation(tRxns); - selO = sparse(1:nT, tRxns, o, nT, nCols2c); dirProb.A = [minProb.A sparse(size(minProb.A,1),nT); ... minProb.c' sparse(1,nT); ... - selO speye(nT)]; + selV sparse(nT,nT) speye(nT)]; dirProb.a = dirProb.A; dirProb.b = [minProb.b; total + max(abs(total)*1e-7, 1e-6); zeros(nT,1)]; dirProb.csense = [minProb.csense 'L' repmat('G',1,nT)]; @@ -658,7 +644,7 @@ %The phase-2c solution is the MIP start: the tight cap on %the total flux leaves the solver's heuristics little room %to find a first solution. - dirProb.start = [res2c.full; max(0, -o.*res2c.full(tRxns))]; + dirProb.start = [res2c.full; max(0, -res2c.full(tRxns))]; res2d = optimizeProb(dirProb,minParams,verbose); [isFeasible2d, isOptimal2d] = checkSolution(res2d); if ~isOptimal2d diff --git a/src/INIT/prepINITModel.m b/src/INIT/prepINITModel.m index d2a27ed3..e129778c 100644 --- a/src/INIT/prepINITModel.m +++ b/src/INIT/prepINITModel.m @@ -194,7 +194,7 @@ %Now, merge all linear dependent rxns - this reduces the size of the model a lot. %Size goes from 11888 rxns to 7922 rxns for human-GEM, so this is an important step -[minModel3,origRxnIds,groupIds,reversedRxns]=mergeLinear(minModel2, {}); +[minModel3,origRxnIds,groupIds]=mergeLinear(minModel2, {}); %we now need to figure out what happened to the essential rxns in the merge step above. %Note that it would be ideal to run checkTasks on the result from mergeLinear, but that doesn't work. @@ -339,7 +339,6 @@ prepData.minModel = scaledMinModel; prepData.refModelWithBM = closeModel(cModel); %do this here so we don't have to do this for each sample prepData.groupIds = groupIds; -prepData.reversedRxns = reversedRxns; %true for a reaction that runs reversed in its merged reaction prepData.essentialMetsForTasks = essentialMetsForTasks; prepData.toIgnoreExch = toIgnoreExch; diff --git a/tests/function_tests/tINIT.m b/tests/function_tests/tINIT.m index 9d22b868..dccb1630 100644 --- a/tests/function_tests/tINIT.m +++ b/tests/function_tests/tINIT.m @@ -481,6 +481,43 @@ function ftINITResolveTiesRanksMergedRxnBySmallestOriginalId(testCase) testCase.verifyFalse(any(ismember({'Rm';'Rc'}, resModel.rxns))); end + function ftINITPinsFluxesInEveryStepButTheLast(testCase) + % ftINIT passes pinFluxes true to ftINITInternalAlg in each step + % whose fluxes a later step reads and false in the last step. A + % wrapper of ftINITInternalAlg earlier on the path records the + % argument. + testCase.assumeMILPSolver(); + global TINIT_INTERNALALG TINIT_PINFLUXES + TINIT_INTERNALALG = @ftINITInternalAlg; + TINIT_PINFLUXES = false(1,0); + d = tempname; mkdir(d); + fid = fopen(fullfile(d,'ftINITInternalAlg.m'),'w'); + fprintf(fid,'%s\n', ... + 'function varargout = ftINITInternalAlg(varargin)', ... + 'global TINIT_INTERNALALG TINIT_PINFLUXES', ... + 'TINIT_PINFLUXES(end+1) = varargin{15};', ... + '[varargout{1:nargout}] = TINIT_INTERNALALG(varargin{:});', ... + 'end'); + fclose(fid); + testCase.addTeardown(@() rmdir(d,'s')); + addpath(d,'-begin'); + testCase.addTeardown(@() rmpath(d)); + testCase.addTeardown(@() clear('global','TINIT_INTERNALALG','TINIT_PINFLUXES')); + + testModel = getTstModelTwoPaths(); + evalc('prepData = prepINITModel(testModel, getTstModelTasks(), {}, false, {}, ''s'');'); + arrayData.genes = testModel.genes; + arrayData.tissues = {'a'}; + arrayData.levels = getExprForRxnScore([-1;-1;-1;-1;5;5]); + arrayData.threshold = 1; + steps = getINITSteps(); + evalc(['resModel = ftINIT(prepData,arrayData.tissues{1},[],[],' ... + 'arrayData,[],steps,true,true,struct(),false,' ... + '''resolveTies'',true);']); + testCase.verifyEqual(TINIT_PINFLUXES, [true(1,numel(steps)-1) false]); + testCase.verifyTrue(all(ismember({'Rn';'Rb'}, resModel.rxns))); + end + function ftINITInternalAlgResolveTiesGivesMinimalFluxes(testCase) % With resolveTies the fluxes come from the least-total-flux % solution for the tie-broken on/off state. The positive @@ -555,7 +592,7 @@ function ftINITInternalAlgResolveTiesWithoutPinFluxesSkipsFluxPhases(testCase) global TINIT_OPTIMIZEPROB_CALLS m = getTstModelTwoSources(); [~,~,~,turnedOn] = ftINITInternalAlg(m,[-1;-1;0;0;0],[],{'E'}, ... - [],false,false,struct(),[],0.1*ones(5,1),false,0.1,true,[],[],false); + [],false,false,struct(),[],0.1*ones(5,1),false,0.1,true,[],false); testCase.verifyEqual(TINIT_OPTIMIZEPROB_CALLS, 3); testCase.verifyEqual(turnedOn(1:2), [true false]); TINIT_OPTIMIZEPROB_CALLS = 0; @@ -564,20 +601,6 @@ function ftINITInternalAlgResolveTiesWithoutPinFluxesSkipsFluxPhases(testCase) testCase.verifyEqual(TINIT_OPTIMIZEPROB_CALLS, 5); end - function ftINITInternalAlgResolveTiesFollowsRxnOrientation(testCase) - % The two-route model of - % ftINITInternalAlgResolveTiesFixesEqualFluxDirections: R1 ranks - % first, so by default R2 carries the backward flux. With - % rxnOrientation -1 for R1, its backward flux in that orientation - % is a forward flux, so R1 runs backward instead. - testCase.assumeMILPSolver(); - [m, scores] = getTstModelTwoRoutes(); - orient = [1;1;-1;1;1;1;1]; - [~,~,~,~,fluxes] = ftINITInternalAlg(m,scores,[],{'E'}, ... - [],false,false,struct(),[],ones(7,1),false,0.1,true,[],orient); - testCase.verifyEqual(fluxes(3:4), [-0.1;0], 'AbsTol', 1e-6); - end - function prepINITModelHandlesModelWithoutReversibleRxns(testCase) % Every reaction is irreversible, so the second simplification % has no reversible reaction to constrain. The merge gives the @@ -589,17 +612,6 @@ function prepINITModelHandlesModelWithoutReversibleRxns(testCase) testCase.verifyEqual(prepData.refModel.rxns, testModel.rxns); end - function prepINITModelReturnsReversedRxnsOfTheMerge(testCase) - % Rz (h[c] => y[c]) and Ra (g[c] => y[c]) both produce y[c], so - % the merge runs Ra reversed in the merged reaction Rz; Rp and - % Rq form a second route between h[c] and g[c]. - testModel = getTstModelFlippedMember(); - evalc('prepData = prepINITModel(testModel, [], {}, false, {}, ''s'');'); - testCase.verifyClass(prepData.reversedRxns, 'logical'); - testCase.verifyEqual(prepData.refModel.rxns(prepData.reversedRxns), {'Ra'}); - testCase.verifyTrue(ismember('Rz', prepData.minModel.rxns)); - end - function ftINITWithTaskRuns(testCase) % A task requiring e[s] from a[s] forces R2 and R7 to be essential. testCase.assumeMILPSolver(); @@ -1122,42 +1134,6 @@ function installOptimizeProbCounter(testCase, fakeTimeLimit) testModel.b = zeros(5,1); end -function testModel = getTstModelFlippedMember() -%h[s] and g[s] are exchanged; h[c] reaches g[c] through y[c] (Rz, Ra) or -%through y2[c] (Rp, Rq). Rz and Ra are both written as producers of y[c]. -testModel = struct(); -testModel.id = 'flippedMember'; -testModel.rxns = {}; -testModel.S = []; -testModel.rev = []; -testModel.mets = {'hs';'hc';'gc';'y1';'y2';'gs'}; -testModel.metNames = {'h';'h';'g';'y';'y2';'g'}; -testModel.comps = {'s';'c'}; -testModel.compNames = testModel.comps; -testModel.metComps = [1;2;2;2;2;1]; -testModel.genes = {'Gz';'Ga';'Gp';'Gq'}; -testModel.grRules = {}; -testModel.rxnGeneMat = []; - -rxnsToAdd = struct(); -rxnsToAdd.rxns = {'EX_h';'Th';'Rz';'Ra';'Rp';'Rq';'Tg';'EX_g'}; -rxnsToAdd.equations = {'<=> h[s]';... - 'h[s] <=> h[c]';... - 'h[c] <=> y[c]';... - 'g[c] <=> y[c]';... - 'h[c] <=> y2[c]';... - 'y2[c] <=> g[c]';... - 'g[c] <=> g[s]';... - 'g[s] <=>'}; -rxnsToAdd.grRules = {'';'';'Gz';'Ga';'Gp';'Gq';'';''}; -evalc('testModel = addRxns(testModel,rxnsToAdd, 3);'); -testModel.c = [0;0;0;0;0;0;0;1]; -testModel.ub = repmat(1000,8,1); -testModel.lb = repmat(-1000,8,1); -testModel.rxnNames = testModel.rxns; -testModel.b = zeros(6,1); -end - function testModel = getTstModelTwoPaths() %a[c] reaches m[c] by Rm, Rc or by Rn, Rb; E, R7 and R8 export e[s]. The %paths are reversible. From b89e453c93fa77f373fb0ac413a70d1c4f260800 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 29 Oct 2026 21:06:02 +0100 Subject: [PATCH 27/27] CHANGELOG: ftINITInternalAlg arguments rankKeys and pinFluxes, phase names of the tie-resolution warning --- CHANGELOG.md | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index d047735b..f197f7a1 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -41,6 +41,11 @@ no step reads its fluxes (new `ftINITInternalAlg` argument `pinFluxes`). The reaction set does not change; the fluxes and `metProduction` of the last step are those of the id-rank phase. +* `ftINITInternalAlg`: new optional positional arguments `rankKeys` (14th, the key each + reaction is sorted by in the id-rank phases, default `model.rxns`) and `pinFluxes` (15th, + default `true`). The "tie resolution did not converge" warning of `resolveTies` also names + `phase2c-minflux` and `phase2d-direction` when the minimal-flux or direction phase does not + prove optimal. * `mergeLinear`: groups are numbered 1, 2, ... in the order of the first reaction of each group in `origRxnIds` (was: the order in which the reactions were merged). * `mergeLinear`: a reaction flipped after the first merge pass is the one marked in