fix: align RAVEN with raven-toolbox (parity review, round 3)

CHANGELOG.md

@@ -1,5 +1,62 @@
 # Changelog
 
+## 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.
+* `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.
+
+### 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
+  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.
+* `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).
+* `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; 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
+  `reversedRxns` (was: the reaction at the same index of the original model).
+
+### 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:

src/INIT/ftINIT.m

@@ -71,10 +71,21 @@
 %     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
-%     gap-filling step (fitTasks' resolveTies) -- see raven-gecko-parity#104
-%     and raven-toolbox's resolve_ties (default false).
+%     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). 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
 %     if set, prove each step to this fixed *absolute* MILP gap in a single
 %     solve, instead of the relative-gap escalation below (stp.MILPParams).
@@ -191,6 +202,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};
@@ -199,6 +222,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;
 
@@ -256,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);
+            [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
@@ -310,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);
+                [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);
 

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,pinFluxes)
 % ftINITInternalAlg
 %	This function runs the MILP for a step in ftINIT.
 %
@@ -44,9 +44,30 @@
 %                   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
-%                   false)
+%                   tie to resolve. The fluxes are then taken 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, 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). 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
+%                   original reaction id of each merge group, so the order
+%                   does not depend on which member names the merged reaction
+%   pinFluxes       if false, resolveTies skips phases 2c and 2d and the
+%                   fluxes and metProduction 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)
@@ -70,9 +91,15 @@
 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
+if nargin < 15 || isempty(pinFluxes)
+    pinFluxes = true;
+end
 if isempty(essentialRxns)
     essentialRxns={};
 end
@@ -453,7 +480,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
@@ -483,6 +513,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 = {};
 
@@ -508,7 +542,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));
@@ -522,6 +556,106 @@
         if isFeasible2b
             %Adopt the tie-resolved solution for the on/off read-out below.
             res = res2b;
+
+            %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
+                %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);
+                        end
+                    end
+                end
+            end
         end
     end
     if ~isempty(unproven)

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
@@ -208,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