OMNI

PURPOSE ^

SYNOPSIS ^

function [OMNISol,bilevelMILPProblem] = OMNI(model, selectedRxnList, options, constrOpt, measOpt, prevSolutions, verbFlag)

DESCRIPTION ^

% ***********************NOT WORKING**************************************

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [OMNISol,bilevelMILPProblem] = OMNI(model, selectedRxnList, options, constrOpt, measOpt, prevSolutions, verbFlag)
0002 %
0003 %% ***********************NOT WORKING**************************************
0004 %
0005 
0006 %function [OMNISol,bilevelMILPProblem] = OMNI(model,selectedRxnList,options,constrOpt,prevSolutions,verbFlag,solutionFileNameTmp)
0007 
0008 %OMNI Run OMNI in the most general form
0009 %
0010 % OMNI(model,selectedRxnList,options,constrOpt,prevSolutions,verbFlag,solutionFileName)
0011 %
0012 %INPUTS
0013 % model                 Structure containing all necessary variables to
0014 %                       describe a stoichiometric model
0015 %   rxns                  Rxns in the model
0016 %   mets                  Metabolites in the model
0017 %   S                     Stoichiometric matrix (sparse)
0018 %   b                     RHS of Sv = b (usually zeros)
0019 %   c                     Objective coefficients
0020 %   lb                    Lower bounds for fluxes
0021 %   ub                    Upper bounds for fluxes
0022 %   rev                    Reversibility of fluxes
0023 % selectedRxnList       List of reactions that can be knocked-out in OMNI
0024 % options               OMNI options
0025 %   numDel                # of bottlenecks
0026 %   numDelSense           Direction of # of bottleneck constraint (G/E/L)
0027 %   vMax                  Max flux
0028 %   solveOMNI             Solve problem within Matlab
0029 %   createGams            Create GAMS input file
0030 %   gamsFile              GAMS input file name
0031 % constrOpt             Explicitly constrained reaction options
0032 %   rxnList               Reaction list
0033 %   values                Values for constrained reactions
0034 %   sense                 Constraint senses for constrained reactions
0035 %                         (G/E/L)
0036 % measOpt               Measured flux options
0037 %   rxnSel                Names of measured reactions
0038 %   values                Flux values of measured reactions
0039 %   weights               Weights for measured fluxes
0040 %
0041 %OPTIONAL INPUTS
0042 % prevSolutions         Previous solutions
0043 % verbFlag              Verbose flag
0044 % solutionFileName      File name for storing temporary solutions
0045 %
0046 %OUTPUTS
0047 % OMNISol               OMNI solution structure
0048 % bilevelMILPProblem    bi-level MILP problem structure used
0049 %
0050 % Markus Herrgard 3/28/05
0051 
0052 % Set these for MILP callbacks
0053 global MILPproblemType;
0054 global selectedRxnIndIrrev;
0055 %global rxnList;
0056 global irrev2rev;
0057 %global solutionFileName;
0058 %global biomassRxnID;
0059 %global OMNIKOrxnList;
0060 %global OMNIObjective;
0061 %global OMNIGrowth;
0062 %global solID;
0063 
0064 if (nargin < 5)
0065     prevSolutions = [];
0066 end
0067 if (nargin < 6)
0068     verbFlag = false;
0069 end
0070 % if (nargin < 7)
0071 %     solutionFileName = 'OMNISolutions.mat';
0072 % else
0073 %     solutionFileName = solutionFileNameTmp;
0074 % end
0075 
0076 % Convert to irreversible rxns
0077 [modelIrrev,matchRev,rev2irrev,irrev2rev] = convertToIrreversible(model);
0078 
0079 % Create the index of the previous KO's suggested by OMNI to avoid obtaining the same
0080 % solution again
0081 selPrevSolIrrev = [];
0082 for i = 1:size(prevSolutions,2)
0083     prevSolRxnList = model.rxns(prevSolutions(:,i)==1);
0084     selPrevSol = ismember(model.rxns,prevSolRxnList);
0085     selPrevSolIrrev(:,i) = selPrevSol(irrev2rev);
0086 end
0087 
0088 [nMets,nRxns] = size(modelIrrev.S);
0089 
0090 % Create matchings for reversible reactions in the set selected for KOs
0091 % This is to ensure that both directions of the reaction are knocked out
0092 selSelectedRxn = ismember(model.rxns,selectedRxnList);
0093 selSelectedRxnIrrev = selSelectedRxn(irrev2rev);
0094 selectedRxnIndIrrev = find(selSelectedRxnIrrev);
0095 cnt = 0;
0096 %prevRxnID = -10;
0097 nSelected = length(selectedRxnIndIrrev);
0098 selRxnCnt = 1;
0099 while selRxnCnt <= nSelected
0100     rxnID = selectedRxnIndIrrev(selRxnCnt);
0101     if (matchRev(rxnID)>0)
0102         cnt = cnt + 1;
0103         selectedRxnMatch(cnt,1) = selRxnCnt;
0104         selectedRxnMatch(cnt,2) = selRxnCnt+1;
0105         selRxnCnt = selRxnCnt + 1;
0106     end
0107     selRxnCnt = selRxnCnt + 1;
0108 end
0109 
0110 % Set inner constraints for the LP
0111 constrOptIrrev = setConstraintsIrrevModel(constrOpt,model,modelIrrev,rev2irrev);
0112 % constrOptIrrev = model;
0113 % constrOptIrrev = [];
0114     
0115 % Set objectives for linear and integer parts
0116 cLinear = zeros(nRxns,1);
0117 cInteger = zeros(sum(selSelectedRxnIrrev),1);
0118 
0119 % Set the correct objective coefficient (not necessary for OMNI)
0120 % targetRxnID = find(ismember(model.rxns,options.targetRxn));
0121 % targetRxnIDirrev = rev2irrev{targetRxnID}(1);
0122 % cLinear(targetRxnIDirrev) = 1;
0123 
0124 % Set measured reaction in objective
0125 sel_meas_rxn = measOpt.rxnSel';
0126 b_meas_rxn = measOpt.values';
0127 wt_meas_rxn = measOpt.weights';
0128 n_m = length(sel_meas_rxn);
0129 
0130 % Create selection vector in the decoupled representation
0131 % This is to ensure that the objective function for measured reversible
0132 % reactions is constructed correctly
0133 sel_m = zeros(nRxns,1);
0134 ord_ir = [];
0135 b_meas_tmp = [];
0136 wt_meas_tmp = [];
0137 for i = 1:n_m
0138     rxn_name = sel_meas_rxn{i};
0139     rxn_id = find(strcmp(model.rxns,rxn_name));
0140     if (~isempty(rxn_id)) % Protect against measured fluxes that are not part of the model
0141         b_meas_tmp = [b_meas_tmp;b_meas_rxn(i)];
0142         wt_meas_tmp = [wt_meas_tmp;wt_meas_rxn(i)];
0143         % Reversible rxns
0144         if (model.rev(rxn_id))
0145             rxn_id_ir = rev2irrev{rxn_id}(1);
0146             sel_m(rxn_id_ir) = 1;
0147             sel_m(rxn_id_ir+1) = -1;
0148         else
0149             % Irrev rxns
0150             rxn_id_ir = rev2irrev{rxn_id};
0151             sel_m(rxn_id_ir) = 1;
0152         end
0153         % Figure out ordering in decoupled representation
0154         ord_ir = [ord_ir rxn_id_ir];
0155     end
0156 end
0157 % Get ordering indices
0158 [tmp,ord_ind] = sort(ord_ir);
0159 % Reorder or create weights
0160 if (sum(wt_meas_rxn) == 0)
0161     measOpts.weights = ones(n_m,1);
0162 else
0163     measOpts.weights = wt_meas_tmp(ord_ind);
0164 end
0165 % Reorder measured flux values
0166 measOpts.values = b_meas_tmp(ord_ind);
0167 
0168 measOpts.rxnSel = sel_m;
0169 
0170 % Create the constraint matrices for the bilevel MILP
0171 bilevelMILPProblem = createBilevelMILPproblem(modelIrrev,cLinear,cInteger,selSelectedRxnIrrev,...
0172     selectedRxnMatch,constrOptIrrev,measOpts,options,selPrevSolIrrev);
0173 
0174 % Initial guess (random)
0175 %bilevelMILPProblem.x0 = round(rand(length(bilevelMILPProblem.c),1));
0176 if isfield(options,'initSolution')
0177     if (length(options.initSolution) > options.numDel | ~all(ismember(options.initSolution,selectedRxnList)))
0178         warning('Initial solution not valid - starting from a random initial solution')
0179         bilevelMILPProblem.x0 = [];
0180     else
0181         % Set initial integer solution
0182         selInitRxn = ismember(model.rxns,options.initSolution);
0183         selInitRxnIrrev = selInitRxn(irrev2rev);
0184         initRxnIndIrrev = find(selInitRxnIrrev);
0185         initIntegerSol = ~ismember(selectedRxnIndIrrev,initRxnIndIrrev);
0186         selInteger = bilevelMILPProblem.vartype == 'B';
0187         [nConstr,nVar] = size(bilevelMILPProblem.A);
0188         bilevelMILPProblem.x0 = nan(nVar,1);
0189         bilevelMILPProblem.x0(selInteger) = initIntegerSol;    
0190         
0191 %         LPproblem.b = bilevelMILPProblem.b - bilevelMILPProblem.A(:,selInteger)*initIntegerSol;
0192 %         LPproblem.A = bilevelMILPProblem.A(:,bilevelMILPProblem.vartype == 'C');
0193 %         LPproblem.c = bilevelMILPProblem.c(bilevelMILPProblem.vartype == 'C');
0194 %         LPproblem.lb = bilevelMILPProblem.lb(bilevelMILPProblem.vartype == 'C');
0195 %         LPproblem.ub = bilevelMILPProblem.ub(bilevelMILPProblem.vartype == 'C');
0196 %         LPproblem.osense = -1;
0197 %         LPproblem.csense = bilevelMILPProblem.csense;
0198 %         LPsol = solveCobraLP(LPproblem);
0199 %
0200 %         bilevelMILPProblem.x0(~selInteger) = LPsol.full;
0201     end
0202 else
0203     bilevelMILPProblem.x0 = [];
0204 end
0205 
0206 % Minimize
0207 bilevelMILPProblem.osense = 1;
0208 
0209 if (verbFlag) 
0210     [nConstr,nVar] = size(bilevelMILPProblem.A);
0211     nInt = length(bilevelMILPProblem.intSolInd);
0212     fprintf('MILP problem with %d constraints %d integer variables and %d continuous variables\n',...
0213         nConstr,nInt,nVar);
0214 end
0215 
0216 bilevelMILPProblem.model = modelIrrev;
0217 
0218 % Set these for CPLEX callbacks
0219 MILPproblemType = 'OMNI';
0220 % rxnList = model.rxns;
0221 % biomassRxnID = find(modelIrrev.c==1);
0222 % solID = 0;
0223 % OMNIObjective = [];
0224 % OMNIGrowth = [];
0225 % OMNIKOrxnList = {};
0226 
0227 % Solve problem
0228 if (options.solveOMNI)
0229     OMNISol = solveCobraMILP(bilevelMILPProblem,'printLevel',0);
0230     if OMNISol.stat~=0
0231         if (~isempty(OMNISol.cont))
0232             OMNISol.fluxes = convertIrrevFluxDistribution(OMNISol.cont(1:length(matchRev)),matchRev);
0233         end
0234         if (~isempty(OMNISol.int))
0235             % Figure out the KO reactions
0236             OMNIRxnInd = selectedRxnIndIrrev(OMNISol.int < 1e-4);
0237             OMNISol.kos = model.rxns(unique(irrev2rev(OMNIRxnInd)));
0238             
0239 %             %sanity check
0240 %             modelTemp = changeRxnBounds(model,OMNISol.kos,0,'b');
0241 %             solTemp = optimizeCbModel(modelTemp);
0242 %             if abs(solTemp.f - OMNISol.obj) > 1e-4
0243 %                 [OMNISol,bilevelMILPProblem] = OMNI(model, selectedRxnList, options, constrOpt, measOpt, prevSolutions, verbFlag);
0244 %                 previous_solutions(:,end+1) = zeros(length(model.rxns),1);
0245 %             end
0246         end
0247     else
0248         OMNISol.fluxes=[];
0249         OMNISol.kos={};
0250     end
0251 else 
0252     OMNISol.rxnList = {};
0253     OMNISol.fluxes = [];
0254 end
0255 
0256 
0257

Generated on Thu 21-Jun-2012 15:39:23 by m2html © 2003