develop changes

src/analysis/exploration/printFluxVector.m

@@ -6,15 +6,15 @@ function printFluxVector(model, fluxData, nonZeroFlag, excFlag, sortCol, fileNam
 %
 % INPUTS:
 %    model:          COBRA model structure
-%    fluxData:       Data matrix/vector (for example, solution.v)
+%    fluxData:       n x k Data matrix/vector (for example, solution.v)
 %
 % OPTIONAL INPUTS:
 %    nonZeroFlag:    Only print nonzero rows (Default = false)
 %    excFlag:        Only print exchange fluxes (Default = false)
 %    sortCol:        Column used for sorting (-1, none; 0, labels; >0, data
 %                    columns;) (Default = -1)
 %    fileName:       Name of output file (Default = [])
-%    headerRow:      Header (Default = [])
+%    headerRow:      k x 1 cell array Header (Default = [])
 %    formulaFlag:    Print reaction formulas (Default = false)
 %    gprFlag:        Print reaction GPR (Default = false)
 %
@@ -47,11 +47,14 @@ function printFluxVector(model, fluxData, nonZeroFlag, excFlag, sortCol, fileNam
     gprFlag = false;
 end
 
+labels = model.rxns;
+labels = pad(labels);
+
 if excFlag
     bool = findExcRxns(model, true, false);
 else
     bool = true(size(model.S,2),1);
-    labels = model.rxns;
+     
 end
 
 if length(nonZeroFlag)==size(model.S,2)
@@ -69,12 +72,11 @@ function printFluxVector(model, fluxData, nonZeroFlag, excFlag, sortCol, fileNam
     if nonZeroFlag
         labels = [model.rxns,model.rxns];
         %only generate the formulas for the nonzero entries
-        formulas = printRxnFormula(model, labels(bool), false, false);
-        labels(bool,2) = formulas;
+        labels(bool,end+1) = printRxnFormula(model, labels(bool), false, false);
     else
-        formulas = printRxnFormula(model, labels, false, false);
-        labels = [labels, formulas];
+        labels(:,end+1) = printRxnFormula(model, labels, false, false);
     end
+    labels(:,end) = pad(labels(:,end));
 end
 
 % Add GPR
@@ -85,6 +87,7 @@ function printFluxVector(model, fluxData, nonZeroFlag, excFlag, sortCol, fileNam
     else
         labels = [labels, model.grRules];
     end
+    labels(:,end) = pad(labels(:,end));
 end
 
 %only print the nonzeros
@@ -93,7 +96,5 @@ function printFluxVector(model, fluxData, nonZeroFlag, excFlag, sortCol, fileNam
     fluxData = fluxData(bool,:);
 end
 
-labels(:,2) = pad(labels(:,2));
-
 %print the labeled data
 printLabeledData(labels, fluxData, 0, sortCol, fileName, headerRow)

src/base/solvers/cplex/buildCplexProblemFromCOBRAStruct.m

@@ -86,7 +86,11 @@
     cplexProblem.Model.ctype = columnVector(Problem.vartype)';
 end
 
-if isfield(Problem,'x0')
-    cplexProblem.Start.x = Problem.x0;
+if isfield(Problem,'basis') || isfield(Problem,'x0')
+    if isfield(Problem,'basis')
+        cplexProblem.Start = Problem.basis;
+    elseif isfield(Problem,'x0')
+        cplexProblem.Start.x = Problem.x0;
+    end
 end
 

src/base/solvers/cplex/setCplexParametersForProblem.m

@@ -1,6 +1,6 @@
 function [cplexProblem,logFile,logToFile] = setCplexParametersForProblem(cplexProblem, cobraParams, solverParams, problemType)
 % Set the parameters for a specific problem from the COBRA Parameter
-% structure and a solver specific parameter structre (latter has
+% structure and a solver specific parameter structure (latter has
 % precedence). The cobra parameters structure contains fields as specified in
 % `getCobraSolverParamsOptionsForType`, while solverParams needs to
 % contain a structure compatible with `setCplexParam`.
@@ -105,6 +105,12 @@
 
 % Set IBM-Cplex-specific parameters. Will overide Cobra solver parameters
 cplexProblem = setCplexParam(cplexProblem, solverParams);
+
+if isfield(cplexProblem,'Start')
+    %https://www.ibm.com/docs/en/icos/12.8.0.0?topic=parameters-advanced-start-switch
+    cplexProblem.Param.advance.Cur = 1;
+end
+
 end
 
 function redirect(outFile,l)

src/base/solvers/msk/parseMskResult.m

@@ -1,4 +1,4 @@
-function [stat,origStat,x,y,z,s,doty] = parseMskResult(res,A,blc,buc,printLevel,param)
+function [stat,origStat,x,y,z,s,doty,bas] = parseMskResult(res,A,blc,buc,printLevel,param)
 %parse the res structure returned from mosek
 
 % initialise variables
@@ -9,6 +9,7 @@
 z = [];
 s = [];
 doty = [];
+bas = [];
 
 if ~exist('printLevel','var')
     printLevel = 0;
@@ -90,6 +91,12 @@
                     % Dual variables to affine conic constraints
                     doty = res.sol.bas.doty;
                 end
+                %https://docs.mosek.com/10.0/toolbox/advanced-hotstart.html
+                bas.skc = res.sol.bas.skc;
+                bas.skx = res.sol.bas.skx;
+                bas.xc = res.sol.bas.xc;
+                bas.xx = res.sol.bas.xx;
+
             case {'PRIMAL_INFEASIBLE_CER','MSK_SOL_STA_PRIM_INFEAS_CER','MSK_SOL_STA_NEAR_PRIM_INFEAS_CER'}
                 stat=0; % infeasible
             case {'DUAL_INFEASIBLE_CER','MSK_SOL_STA_DUAL_INFEAS_CER','MSK_SOL_STA_NEAR_DUAL_INFEAS_CER'}

src/base/solvers/solveCobraLP.m

@@ -705,25 +705,41 @@
         %              blx and bux. Note -inf is allowed in blc and blx.
         %              Similarly, inf is allowed in buc and bux.
 
-        if isempty(csense)
-            % assumes all equality constraints
-            % [res] = msklpopt(c,a,blc,buc,blx,bux,param,cmd)
-            [res] = msklpopt(osense * c, A, b, b, lb, ub, param, cmd);
-        else
-            blc = b;
-            buc = b;
-            buc(csense == 'G') = inf;
-            blc(csense == 'L') = -inf;
-            % [res] = msklpopt(       c,a,blc,buc,blx,bux,param,cmd)
+
+        blc = b;
+        buc = b;
+        buc(csense == 'G') = inf;
+        blc(csense == 'L') = -inf;
+
+        if 0
             [res] = msklpopt(osense * c, A, blc, buc, lb, ub, param, cmd);
-%             res.sol.itr
-%             min(buc(csense == 'E')-A((csense == 'E'),:)*res.sol.itr.xx)
-%             min(A((csense == 'E'),:)*res.sol.itr.xx-blc(csense == 'E'))
-%             pasue(eps)
+            %             res.sol.itr
+            %             min(buc(csense == 'E')-A((csense == 'E'),:)*res.sol.itr.xx)
+            %             min(A((csense == 'E'),:)*res.sol.itr.xx-blc(csense == 'E'))
+            %             pasue(eps)
+
+        else
+            prob.c = osense * c;
+            prob.a = A;
+            prob.blc     = blc;
+            prob.buc     = buc;
+            prob.blx     = lb;
+            prob.bux     = ub;
+
+            if ~isempty(basis)
+                prob.sol.bas.skc  = basis.skc;
+                prob.sol.bas.skx  = basis.skx;
+                prob.sol.bas.xc   = basis.xc;
+                prob.sol.bas.xx   = basis.xx;
+            end
+
+            % Use the primal simplex optimizer.
+            param.MSK_IPAR_OPTIMIZER = 'MSK_OPTIMIZER_PRIMAL_SIMPLEX';
+            [rcode,res] = mosekopt(cmd,prob,param);
         end
-
+                
         %parse mosek result structure
-        [stat,origStat,x,y,w,s,~] = parseMskResult(res,A,blc,buc,problemTypeParams.printLevel,param);
+        [stat,origStat,x,y,w,s,~,basis] = parseMskResult(res,A,blc,buc,problemTypeParams.printLevel,param);
         if stat ==1
             f=c'*x;
         else
@@ -1276,6 +1292,7 @@
         [CplexLPproblem, logFile, logToFile] = setCplexParametersForProblem(CplexLPproblem,problemTypeParams,solverParams,'LP');
         %logToFile=0;
 
+
         % optimize the problem
         CplexLPproblem.solve();
 
@@ -1299,6 +1316,10 @@
             y = osense*CplexLPproblem.Solution.dual;
             %res1 = A*solution.full + solution.slack - b;
             s = b - A * x; % output the slack variables
+
+            %save basis also
+            basis = CplexLPproblem.Start;
+
         elseif origStat == 2 ||   origStat == 20
             stat = 2; %unbounded
         elseif origStat == 3
@@ -1374,6 +1395,8 @@
         if exist([pwd filesep 'clone2_' labindex '.log'],'file')
             delete([pwd filesep 'clone2_' labindex '.log'])
         end
+
+
     case 'lindo'
         %%
         error('The lindo interface is obsolete.');

src/dataIntegration/XomicsToModel/XomicsToModel.m

@@ -31,15 +31,19 @@
 %    * .C - `k x n` Left hand side of C*v <= d
 %    * .d - `k x 1` Right hand side of C*v <= d
 %    * .dsense - `k x 1` character array with entries in {L,E,G}
+%    * .beta - A scalar weight on minimisation of one-norm of internal fluxes. Default 1e-4. 
+%              Larger values increase the incentive to find a flux vector to be thermodynamically feasibile in thermoKernel and decrease the incentive
+%              to search the steady state solution space for a flux vector that results in certain reactions and metabolites to be active and present,
+%              respectively.
 %
 %%  specificData:  A structure containing context-specific data:
 %
 %   * .activeGenes - cell array of Entrez ID of genes that are known to be active based on the bibliomic data (Default: empty).
 %   * .inactiveGenes - cell array of Entrez ID of genes known to be inactive based on the bibliomics data (Default: empty).
 %
-%   * .activeReactions -cell array of reactions know to be active based on bibliomic data (Default: empty).
-%   * .inactiveReactions - cell array of reactions know to be inactive based on bibliomic data (Default: empty).
-
+%   * .activeReactions -cell array of reaction identifiers know to be active based on bibliomic data (Default: empty).
+%   * .inactiveReactions - cell array of reaction identifiers know to be inactive based on bibliomic data (Default: empty).
+i
 %   * .coupledRxns -Table containing information about the coupled reactions. This includes the coupled reaction identifier, the
 %                    list of coupled reactions, the coefficients of those reactions, the constraint, the sense or the directionality of the constraint,
 %                    and the reference (Default: empty).
@@ -556,7 +560,7 @@
 end
 % Fix demand reaction names
 %DM_CE5026[c] -> Demand for 5-S-Glutathionyl-L-Dopa
-rxnNamesToFixList = model.rxns(contains(model.rxns, {'DM_', 'EX_', 'sink_'}) & strcmp(model.rxns, model.rxnNames));
+rxnNamesToFixList = model.rxns(contains(model.rxns, {'DM_', 'EX_', '_'}) & strcmp(model.rxns, model.rxnNames));
 for i = 1:length(rxnNamesToFixList)
     rxnBool = strcmp(model.rxns, rxnNamesToFixList{i});
     metBool = model.S(:, rxnBool) ~= 0;
@@ -1120,7 +1124,7 @@
 %% 9. Close sink and demand reactions - Set non-core sinks and demands to inactive
 coreRxnBool = ismember(model.rxns, coreRxnAbbr);
 
-model.lbpreSinkDemandOff = model.lb;
+model.model = model.lb;
 model.ubpreSinkDemandOff = model.ub;
 
 % Identify reaction type

src/dataIntegration/XomicsToModel/thermoKernel/thermoKernel.m

@@ -16,13 +16,16 @@
 %
 % OPTIONAL INPUTS:
 %    model:             (optional fields)
-%          * b - `m x 1` change in concentration with time
-%          * csense - `m x 1` character array with entries in {L,E,G}
-%          * osenseStr: Maximize ('max')/minimize ('min') (opt, default = 'max') linear part of the objective.
-%          * C - `k x n` Left hand side of C*v <= d
-%          * d - `k x n` Right hand side of C*v <= d
-%          * dsense - `k x 1` character array with entries in {L,E,G}
-%          * beta - scalar  trade-off parameter on minimisation of one-norm of internal fluxes. Increase to incentivise thermodynamic feasibility in optCardThermo
+%          * .b - `m x 1` change in concentration with time
+%          * .csense - `m x 1` character array with entries in {L,E,G}
+%          * .osenseStr: Maximize ('max')/minimize ('min') (opt, default = 'max') linear part of the objective.
+%          * .C - `k x n` Left hand side of C*v <= d
+%          * .d - `k x n` Right hand side of C*v <= d
+%          * .dsense - `k x 1` character array with entries in {L,E,G}
+%          * .beta - A scalar weight on minimisation of one-norm of internal fluxes. Default 1e-4. 
+%                    Larger values increase the incentive to find a flux vector to be thermodynamically feasibile in each iteration of optCardThermo 
+%                    and decrease the incentive to search the steady state solution space for a flux vector that results in certain reactions and
+%                    metabolites to be active and present, respectively.
 %
 %    activeInactiveRxn: - `n x 1` with entries {1,-1, 0} depending on whether a reaction must be active, inactive, or unspecified respectively.
 %    rxnWeights:        - `n x 1` real valued penalties on zero norm of reaction flux, negative to promote a reaction to be active, positive 
@@ -42,9 +45,17 @@
 %                   * .normalizeZeroNormWeights - {(0),1}, normalises zero norm weights
 %                                                 rxnWeights  = rxnWeights./sum(abs(rxnWeights));
 %                                                 metWeights  = metWeights./sum(abs(metWeights));
-%                   * .param.removeOrphanGenes - {(1),0}, removes orphan genes from thermoModel
-%                   * .formulation - mathematical formulation of thermoKernel algorithm (use default unless expert user)
-% 
+%                   * .rxnWeightsConsistentWithMetWeights {(0),1} If true and metWeights are provided, make the corresponding reaction weights consistent
+%                   * .metWeightsConsistentWithRxnWeights {(0),1} If true and rxnWeights are provided, make the corresponding metabolite weights consistent
+%                   * .acceptRepairedFlux {(1),0} If true, a post processing step after each inner iteration minimises the absolute value of internal reaction flux, 
+%                                                 while (a) all exchange fluxes are kept constant, and (b) no internal flux is allowed to change direction or increase in size.   
+%                   * .relaxBounds {(0),1} If true, allow internal bounds forcing non-zero flux to be relaxed as minimising absolute value of internal 
+%                                          fluxes is only guarunteed to return a thermodynamically feasible flux if such bounds can be relaxed.
+%                   * .removeOrphanGenes - {(1),0}, removes orphan genes from thermoModel
+%                   * .formulation - mathematical formulation of thermoKernel algorithm (Default is 'pqzwrs'. Do not change unless expert user.)
+%                   * .plotThermoKernelStats {(0),1} generates a figure with confusion matrices comparing anticipated vs actual metabolites and reactions in the extracted model
+%                   * .plotThermoKernelWeights {(0),1} generates a figure displaying the weights given to actual and anticipated but omitted metabolites and reactions in the extracted model
+
 % OUTPUTS:
 %   thermoModel:           thermodynamically consistent model extracted from input model
 %   thermoModelMetBool:   `m` x 1 boolean vector of thermodynamically consistent `mets` in input model