Skip to content

Commit

Permalink
Merge pull request #2097 from opencobra/develop
Browse files Browse the repository at this point in the history 
Develop
  • Loading branch information
@rmtfleming
rmtfleming authored Dec 15, 2022
2 parents 81d019b + eee472d commit 44a3e1e
Show file tree
Hide file tree
Showing 4 changed files with 54 additions and 68 deletions.
2 changes: 1 addition & 1 deletion ...tiSpecies/microbiomeModelingToolbox/pairwiseInteractionModeling/computeParetoOptimality.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -71,7 +71,7 @@
model.osenseStr = 'max';
[solution] = solveCobraLP(buildLPproblemFromModel(model));
dmaxRxn1 = solution.obj;
model = changeObjective(model, rxn2);
model = changeObjective(model, rxn1);
model.osenseStr = 'min';
[solution] = solveCobraLP(buildLPproblemFromModel(model));
dminRxn1 = solution.obj;
Expand Down
107 changes: 45 additions & 62 deletions src/analysis/thermo/findThermoConsistentFluxSubset/findThermoConsistentFluxSubset.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,7 @@
%
% param.formulation: mathematical formulation of inner iteration
% param.relaxBounds: Relax bounds that don't include zero. Default is false.
% param.nMax
%
% OUTPUT
% thermoFluxConsistentMetBool m x 1 boolean vector indicating thermodynamically flux consistent mets
Expand Down Expand Up @@ -232,79 +233,61 @@
%recalculate size of S
[nMet,nRxn]=size(model.S);

%secondary removal is on by default
if param.secondaryRemoval
% Finds the subset of `S` that is stoichiometrically consistent using
% an iterative cardinality optimisation approach
%heuristically identify exchange reactions and metabolites exclusively
%involved in exchange reactions
if ~isfield(model,'SIntRxnBool') || ~isfield(model,'SIntMetBool')
if isfield(model,'mets')
%attempts to finds the reactions in the model which export/import from the model
%boundary i.e. mass unbalanced reactions
%e.g. Exchange reactions
% Demand reactions
% Sink reactions
model = findSExRxnInd(model,[],param.printLevel-1);
else
model.SIntMetBool=true(size(model.S,1),1);
model.SIntRxnBool=true(size(model.S,2),1);
end
else
if length(model.SIntMetBool)~=size(model.S,1) || length(model.SIntRxnBool)~=size(model.S,2)
model = findSExRxnInd(model,[],param.printLevel-1);
end
end

%extract the stoichiometrically consistent subset
if param.secondaryRemoval && (~isfield(model,'SConsistentRxnBool') || ~isfield(model,'SConsistentMetBool') || length(model.SConsistentMetBool)~=size(model.S,1) || length(model.SConsistentRxnBool)~=size(model.S,2))

SIntRxnBool = model.SIntRxnBool;

massBalanceCheck=0;
[SConsistentMetBool, SConsistentRxnBool, ~, ~, ~, ~, model] = findStoichConsistentSubset(model,massBalanceCheck);
[SConsistentMetBool, SConsistentRxnBool, SInConsistentMetBool, SInConsistentRxnBool, unknownSConsistencyMetBool, unknownSConsistencyRxnBool, ~, model]...
= findStoichConsistentSubset(model, massBalanceCheck, param.printLevel-1);

if param.printLevel>0
fprintf('%10u%10u%s\n',nnz(~SConsistentMetBool),nnz(model.SIntRxnBool & ~SConsistentRxnBool),' metabolites and reactions removed by findStoichConsistentSubset called from findThermoConsistentFluxSubset.');
fprintf('%10u%%s\n',nnz(~SConsistentMetBool), ' metabolites removed by findStoichConsistentSubset called from findThermoConsistentFluxSubset.');
end
if param.printLevel>0
fprintf('%10u%%s\n',nnz(SIntRxnBool & ~SConsistentRxnBool),' reactions removed by findStoichConsistentSubset called from findThermoConsistentFluxSubset.');
end

end

%extract the flux consistent subset
if param.secondaryRemoval && (~isfield(model,'fluxConsistentRxnBool') || ~isfield(model,'fluxConsistentMetBool') || length(model.fluxConsistentMetBool)~=size(model.S,1) || length(model.fluxConsistentRxnBool)~=size(model.S,2))
if isfield(model,'C')
fluxConsistentParam.method='fastcc';%can handle additional constraints
else
fluxConsistentParam.method='null_fastcc';
end
[fluxConsistentMetBool, fluxConsistentRxnBool,~,~,model]= findFluxConsistentSubset(model,fluxConsistentParam);
if param.printLevel>0
fprintf('%10u%10u%s\n',nnz(~fluxConsistentMetBool),nnz(~fluxConsistentRxnBool),' metabolites and reactions removed by findFluxConsistentSubset called from findThermoConsistentFluxSubset.');
end

% if 1
% metRemoveList1 = model.mets(~fluxConsistentMetBool);
%
% [model, metRemoveList2] = removeRxns(model, model.rxns(~fluxConsistentRxnBool),'metRemoveMethod','exclusive','ctrsRemoveMethod','inclusive');
% if length(intersect(metRemoveList1,metRemoveList2))~=length(metRemoveList1)
% error('inconsistent metabolite removal')
% end
%
% else
% model.S = model.S(model.fluxConsistentMetBool,model.fluxConsistentRxnBool);
% model.lb = model.lb(model.fluxConsistentRxnBool);
% model.ub = model.ub(model.fluxConsistentRxnBool);
% model.c = model.c(model.fluxConsistentRxnBool);
% model.SConsistentRxnBool = model.SConsistentRxnBool(model.fluxConsistentRxnBool);
% model.SConsistentMetBool = model.SConsistentMetBool(model.fluxConsistentMetBool);
%
% [nMet,nRxn]=size(model.S);
%
% if isfield(model,'b')
% model.b = model.b(model.fluxConsistentMetBool);
% else
% model.b = zeros(nMet,1);
% end
% % assume constraint S*v = b if csense not provided
% if isfield(model, 'csense')
% model.csense = model.csense(model.fluxConsistentMetBool);
% else
% % if csense is not declared in the model, assume that all
% % constraints are equalities.
% model.csense(1:nMet, 1) = 'E';
% end
%
% if isfield(model,'C')
% %constraints inlusively involved in flux inconsistent reactions are deemed flux inconsistent also
% fluxConsistentConstraintBool = getCorrespondingRows(model.C,true(size(model.C,1),1),model.fluxConsistentRxnBool,'inclusive');
%
% if ~all(fluxConsistentConstraintBool)
% fprintf('%s\n',[num2str(nnz(~fluxConsistentConstraintBool)) ' model.C constraints removed'])
%
% model.C = model.C(fluxConsistentConstraintBool,model.fluxConsistentRxnBool);
% model.d = model.d(fluxConsistentConstraintBool,1);
% model.dsense = model.dsense(fluxConsistentConstraintBool,1);
% if isfield(model,'ctrs')
% model.ctrs = model.ctrs(fluxConsistentConstraintBool,1);
% end
% if isfield(model,'ctrNames')
% model.ctrNames = model.ctrNames(fluxConsistentConstraintBool,1);
% end
% end
% end
% end
fluxConsistentParam.printLevel=param.printLevel-1;
[fluxConsistentMetBool, fluxConsistentRxnBool, ~, ~,~,model] = findFluxConsistentSubset(model,fluxConsistentParam);

if any(~fluxConsistentRxnBool)
fprintf('%u%s\n',nnz(~fluxConsistentRxnBool), ' flux inconsistent reactions removed by findFluxConsistentSubset called from findThermoConsistentFluxSubset.')
end

if any(~fluxConsistentMetBool)
fprintf('%u%s\n',nnz(~fluxConsistentMetBool), ' flux inconsistent metabolites removed by findFluxConsistentSubset called from findThermoConsistentFluxSubset.')
end
end

%%
Expand Down
6 changes: 2 additions & 4 deletions src/base/solvers/solveCobraQP.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1038,8 +1038,7 @@

% evaluate the optimality condition 1
if tmp1 > problemTypeParams.feasTol * 1e2
disp(solution.origStat)
warning(['[' solver '] Primal optimality condition in solveCobraQP not satisfied, residual = ' num2str(tmp1) ', while feasTol = ' num2str(problemTypeParams.feasTol)])
fprintf('%s\n',['[' solver '] reports ' solution.origStat ' but Primal optimality condition in solveCobraQP not satisfied, residual = ' num2str(tmp1) ', while feasTol = ' num2str(problemTypeParams.feasTol)])
else
if problemTypeParams.printLevel > 0
fprintf(['\n > [' solver '] Primal optimality condition in solveCobraQP satisfied.']);
Expand All @@ -1053,8 +1052,7 @@

% evaluate the optimality condition 2
if tmp2 > problemTypeParams.optTol * 1e2
disp(solution.origStat)
warning(['[' solver '] Dual optimality condition in solveCobraQP not satisfied, residual = ' num2str(tmp2) ', while optTol = ' num2str(problemTypeParams.optTol)])
fprintf('%s\n',['[' solver '] reports ' solution.origStat ' but Dual optimality condition in solveCobraQP not satisfied, residual = ' num2str(tmp2) ', while optTol = ' num2str(problemTypeParams.optTol)])
else
if problemTypeParams.printLevel > 0
fprintf(['\n > [' solver '] Dual optimality condition in solveCobraQP satisfied.\n']);
Expand Down
7 changes: 6 additions & 1 deletion src/base/utilities/mapAontoB.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -58,7 +58,12 @@
end
Bout=Bout(ileft~=0,:);
else
Bout(LIBkey,:) = Ain(LOCAkey(LOCAkey~=0),:);
Bout = Bin;
if size(Bout,2)==size(Ain,2)
Bout(LIBkey,:) = Ain(LOCAkey(LOCAkey~=0),:);
else
error('size(Bout,2)~=size(Ain,2)')
end
end
else
switch classAin
Expand Down

0 comments on commit 44a3e1e

Please sign in to comment.