optimise cardinality reduced regularisation

Author: rmtfleming

binary

@@ -35,6 +35,7 @@ function plotConfMat(varargin)
 % plotting the colors
 imagesc(confpercent);
 title(sprintf('Accuracy: %.2f%%', 100*trace(confmat)/sum(confmat(:))));
+%ylabel('Predicted Class'); xlabel('Target Class');
 ylabel('Predicted Class'); xlabel('Target Class');
 
 % set the colormap

external/analysis/plotConfMat/plotConfMat.m

@@ -224,7 +224,7 @@ function initCobraToolbox(updateToolbox)
         %This means that the checked-out commit -- which is the one that the super-project (core) needs -- is not associated with a local branch name.
         %[status_gitSubmodule, result_gitSubmodule] = system(['git submodule update --init --remote --no-fetch ' depthFlag]);%old
         %[status_gitSubmodule, result_gitSubmodule] = system(['git submodule foreach git submodule update --init --recursive']);% 23/9/21 RF submodules point to master
-        [status_gitSubmodule, result_gitSubmodule] = system('git submodule update --init --recursive');% 23/9/21 RF submodules point to master, don't pull in remote changes
+        [status_gitSubmodule, result_gitSubmodule] = system('git submodule update --init --recursive --depth 1');% 23/9/21 RF submodules point to master, don't pull in remote changes
         %[status_gitSubmodule, result_gitSubmodule] = system('git submodule foreach git checkout master');
         [status_gitSubmodule, result_gitSubmodule] = system('git submodule foreach git checkout master');% 30/9/21 RF submodules point to master, don't pull in remote changes
 

initCobraToolbox.m

@@ -45,16 +45,16 @@
 %                         * ctrs `k x 1` Cell Array of Strings giving IDs of the coupling constraints
 %
 %                         * dsense - `k x 1` character array with entries in {L,E,G}
-%                         * g0 - `n x 1` weights on zero norm
-%                         * g1 - `n x 1` weights on one norm
+%                         * g0 - `n x 1` weights on zero norm, where positive is minimisation, negative is maximisation, zero is neither.
+%                         * g1 - `n x 1` weights on one norm, where positive is minimisation, negative is maximisation, zero is neither.
 %                         * g2 - `n x 1` weights on two norm
 %
 %    osenseStr:         Maximize ('max')/minimize ('min') (opt, default =
 %                       'max') linear part of the objective. Nonlinear
 %                       parts of the objective are always assumed to be
 %                       minimised.
 %
-%    minNorm:           {(0), 'one', 'zero', > 0 , n x 1 vector}, where `[m,n]=size(S)`;
+%    minNorm:           {(0), 'one', 'zero', > 0 , n x 1 vector, 'optimizeCardinality'}, where `[m,n]=size(S)`;
 %                       0 - Default, normal LP
 %                       'one'  Minimise the Taxicab Norm using LP.
 %
@@ -66,6 +66,7 @@
 %                               ~& lb \leq v \leq ub
 %
 %                       A LP solver is required.
+
 %                       'zero' Minimize the cardinality (zero-norm) of v
 %
 %                       .. math::
@@ -86,6 +87,9 @@
 %                          http://dx.doi.org/10.1016/j.ejor.2014.11.031
 %                          A LP solver is required.
 %
+%                       'optimizeCardinality' as for 'zero' option but uses
+%                       model.g0 - `n x 1` weights on zero norm, where positive is minimisation, negative is maximisation, zero is neither.
+%
 %                       The remaining options work only with a valid QP solver:
 %
 %                       > 0    Minimises the squared Euclidean Norm of internal fluxes.
@@ -107,6 +111,8 @@
 %                          s.t. ~& S v = b \\
 %                               ~& c^T v = f \\
 %                               ~& lb \leq v \leq ub
+%
+
 %
 %    allowLoops:        {0,(1)} If false, then instead of a conventional FBA,
 %                       the solver will run an MILP version which does not allow

src/analysis/FBA/optimizeCbModel.m

@@ -76,6 +76,7 @@
 stop = false;
 solution.stat = 1;
 
+feasTol = getCobraSolverParams('LP', 'feasTol');
 
 if exist('param','var')
     if isfield(param,'excludedReactions') == 0
@@ -95,11 +96,6 @@
         param.nbMaxIteration = 1000;
     end
 
-    if isfield(param,'epsilon') == 0
-        feasTol = getCobraSolverParams('LP', 'feasTol');
-        param.epsilon = feasTol*100;
-    end
-
     if isfield(param,'gamma0') == 0
         param.gamma0 = 0;    %trade-off parameter of l0 part v
     end
@@ -108,27 +104,29 @@
         %Small positive penalty seems to be important to ensure algorithm
         %is numerically stable for small networks.
         %TODO: why?
-        param.gamma1 = 1e-6;     %trade-off parameter of l1 part v
+        param.gamma1 = feasTol*100;     %trade-off parameter of l1 part v
+        
     end
 
     if isfield(param,'lambda0') == 0
-        param.lambda0 = 10;   %trade-off parameter of l0 part of r
+        param.lambda0 = 1;   %trade-off parameter of l0 part of r
     end
 
     if isfield(param,'lambda1') == 0
-        param.lambda1 = 1;    %trade-off parameter of l1 part of r
+        param.lambda1 = feasTol*100;    %trade-off parameter of l1 part of r
     end
 
     if isfield(param,'alpha0') == 0
-        param.alpha0 = 10;   %trade-off parameter of l0 part of p and q
+        param.alpha0 = 1;   %trade-off parameter of l0 part of p and q
     end
 
     if isfield(param,'alpha1') == 0
-        param.alpha1 = 1;    %trade-off parameter of l1 part of p and q
+        param.alpha1 = feasTol*100;    %trade-off parameter of l1 part of p and q
+   
     end
-
-    if isfield(param,'epsilon') == 0
-        param.epsilon = 10e-6; %stopping criterion
+  
+    if isfield(param,'epsilon') == 0  %stopping criterion
+        param.epsilon = feasTol*100;
     end
 
     if isfield(param,'theta') == 0
@@ -137,7 +135,13 @@
 end
 
 if 0
-    param
+    param.gamma1 = 0; % caffeine
+    param.lambda1 = 0; % caffeine
+    param.alpha1 = 0; % caffeine
+end
+
+if param.printLevel>1
+    disp(param)
 end
 
 [nbMaxIteration,epsilon,theta]      = deal(param.nbMaxIteration,param.epsilon,param.theta);
@@ -208,7 +212,7 @@
         v_bar = -gamma1*sign(v) + sign(v)*gamma0*theta;
 
         r(abs(r) < one_over_theta) = 0;
-        r_bar = -lamda1*sign(r) + sign(r)*lambda0*theta;
+        r_bar = -lambda1*sign(r) + sign(r)*lambda0*theta;
 
         p(p < one_over_theta) = 0;
         p_bar = alpha1*sign(p) + sign(p)*alpha0*theta;
@@ -253,7 +257,7 @@
                 if nbIteration==0
                     fprintf('%5s%12s%12s%12s%12s%10s%10s%10s%10s\n','itn','obj','obj_old','err(obj)','err(x)','card(v)','card(r)','card(p)','card(q)');
                 end
-                fprintf('%5u%12.5g%12.5g%12.5g%12.5g%10u%10u%10u%10u\n',nbIteration,obj_new,obj_old,error_obj,error_x,nnz(v),nnz(r),nnz(p),nnz(q));
+                fprintf('%5u%12.5g%12.5g%12.5g%12.5g%10u%10u%10u%10u\n',nbIteration,obj_new,obj_old,error_obj,error_x,nnz(abs(v)>feasTol),nnz(abs(r)>feasTol),nnz(p>feasTol),nnz(q>feasTol));
             end
 
             if (error_x < epsilon) || (error_obj < epsilon)
@@ -283,7 +287,7 @@
                 if nbIteration==0
                     fprintf('%5s%12s%12s%12s%12s%10s%10s%10s%10s\n','itn','obj','obj_old','err(obj)','err(x)','card(v)','card(r)','card(p)','card(q)');
                 end
-                fprintf('%5u%12.5g%12.5g%12.5g%12.5g%10u%10u%10u%10u\n',nbIteration,obj_new,obj_old,error_obj,error_x,nnz(v),nnz(r),nnz(p),nnz(q));
+                fprintf('%5u%12.5g%12.5g%12.5g%12.5g%10u%10u%10u%10u\n',nbIteration,obj_new,obj_old,error_obj,error_x,nnz(abs(v)>feasTol),nnz(abs(r)>feasTol),nnz(p>feasTol),nnz(q>feasTol));
             end
             stop = true;
     end

src/analysis/relaxedFBA/relaxFBA_cappedL1.m

@@ -39,48 +39,40 @@
 %    param:    Structure optionally containing the relaxation parameters:
 %
 %                      * .internalRelax:
-%
 %                        * 0 = do not allow to relax bounds on internal reactions
 %                        * 1 = do not allow to relax bounds on internal reactions with finite bounds
 %                        * {2} = allow to relax bounds on all internal reactions
 %
 %                      * .exchangeRelax:
-%
 %                        * 0 = do not allow to relax bounds on exchange reactions
 %                        * 1 = do not allow to relax bounds on exchange reactions of the type [0,0]
 %                        * {2} = allow to relax bounds on all exchange reactions
 %
 %                      * .steadyStateRelax:
-%
 %                        *    0 = do not allow to relax the steady state constraint S*v = b
 %                        *  {1} = allow to relax the steady state constraint S*v = b
 %
 %                      * .toBeUnblockedReactions - nRxns x 1 vector indicating the reactions to be unblocked
-%
 %                        * toBeUnblockedReactions(i) = 1 : impose v(i) to be positive
 %                        * toBeUnblockedReactions(i) = -1 : impose v(i) to be negative
 %                        * toBeUnblockedReactions(i) = 0 : do not add any constraint (default)
 %
 %                      * .excludedReactions - nRxns x 1 bool vector indicating the reactions to be excluded from relaxation
-%
 %                        * excludedReactions(i) = false : allow to relax bounds on reaction i (default)
 %                        * excludedReactions(i) = true : do not allow to relax bounds on reaction i 
 %
 %                      * .excludedReactionLB - nRxns x 1 bool vector indicating
 %                      the reactions with lower bounds to be excluded from
 %                      relaxation (overridden by excludedReactions)
-%
 %                        * excludedReactionLB(i) = false : allow to relax lower bounds on reaction i (default)
 %                        * excludedReactionLB(i) = true : do not allow to relax lower bounds on reaction i 
 %
 %                      * .excludedReactionUB - nRxns x 1 bool vector indicating
 %                      the reactions with upper bounds to be excluded from relaxation (overridden by excludedReactions)
-%
 %                        * excludedReactionUB(i) = false : allow to relax upper bounds on reaction i (default)
 %                        * excludedReactionUB(i) = true : do not allow to relax upper bounds on reaction i 
 %
 %                      * .excludedMetabolites - nMets x 1 bool vector indicating the metabolites to be excluded from relaxation
-%
 %                        * excludedMetabolites(i) = false : allow to relax steady state constraint on metabolite i (default)
 %                        * excludedMetabolites(i) = true : do not allow to relax steady state constraint on metabolite i
 %
@@ -99,9 +91,8 @@
 %                     the algorithm is useful when trying different values
 %                     of theta to start with the appropriate parameter
 %                     giving the lowest cardinality solution.
-%                     * .relaxedPrintLevel (Default = 0) Printing information on relaxed reactions
-%                     * .maxRelaxR (Default = 1e4), maximum relaxation
-%                     of any bound or equality constraint permitted
+%                     * .relaxedPrintLevel (Default = 0) Printing information on relaxed reaction bounds and steady state constraints
+%                     * .maxRelaxR (Default = 1e4), maximum relaxation of any bound or equality constraint permitted
 %
 % OUTPUT:
 %    solution:       Structure containing the following fields:
@@ -127,6 +118,8 @@
     issueConfirmationWarning('relaxedFBA ignores additional variables defined in the model (model field .E)!')
 end
 
+feasTol = getCobraSolverParams('LP', 'feasTol');
+
 if ~isfield(model,'S')
     %relax generic LP problem ,e.g.
 %          A: [1663×2942 double]
@@ -159,7 +152,7 @@
     param.minLB = min(-max(model.ub),min(model.lb));
 end
 if ~isfield(param,'maxRelaxR')
-    param.maxRelaxR = 1e4; %TODO - check this for multiscale models
+    param.maxRelaxR = 1/feasTol; %TODO - check this for multiscale models
 end
 if ~isfield(param,'printLevel')
     param.printLevel = 0; %TODO - check this for multiscale models
@@ -217,7 +210,7 @@
 end
 %TODO properly incorporate inf bounds
 %add a large finite lower bound here
-model.lb(model.lb==-inf) = -1e6;
+model.lb(model.lb==-inf) = -1/feasTol;
 %use this to override some other assignment
 excludedReactionLBTmp=param.excludedReactionLB | model.lb==-inf;
 
@@ -228,7 +221,7 @@
     param.excludedReactionUB = false(nRxns,1);
 end
 %add a large finite upper bound here
-model.ub(model.ub==inf) = 1e6;
+model.ub(model.ub==inf) = 1/feasTol;
 %use this to override some other assignment
 excludedReactionUBTmp=param.excludedReactionUB | model.ub==inf;
 
@@ -250,8 +243,7 @@
 end
 
 if isfield(param,'epsilon') == 0
-    feasTol = getCobraSolverParams('LP', 'feasTol');
-    param.epsilon=feasTol*100;
+    param.epsilon=feasTol;
 end
 
 if isfield(param,'theta') == 0
@@ -353,10 +345,10 @@
 
 %set global parameters on zero norm if they do not exist
 if ~isfield(param,'alpha') && ~isfield(param,'alpha0')
-    param.alpha = 10; %weight on relaxation of bounds of reactions
+    param.alpha = 1; %weight on relaxation of bounds of reactions
 end
 if ~isfield(param,'lambda') && ~isfield(param,'lambda0')
-    param.lambda = 10;  %weight on relaxation of steady state constraints
+    param.lambda = 1;  %weight on relaxation of steady state constraints
 end
 if ~isfield(param,'gamma') && ~isfield(param,'gamma0')
     %default should not be to aim for zero norm flux vector if the problem is infeasible at the begining 
@@ -376,15 +368,14 @@
 
 %set local paramters on one norm for capped L1
 if ~isfield(param,'alpha1')
-    param.alpha1 = param.alpha0/10; 
+    param.alpha1 = feasTol; 
 end
 if ~isfield(param,'lambda1')
-    param.lambda1 = param.lambda0/10;    
+    param.lambda1 = feasTol;    
 end
 if ~isfield(param,'gamma1')
-    %always include some regularisation on the flux rates to keep it well
-    %behaved
-    param.gamma1 = 1e-6 + param.gamma0/10;   
+    %some regularisation on the flux rates to keep it well behaved
+    param.gamma1 = feasTol;   
 end
 
 %Combine excludedReactions with internalRelax and exchangeRelax
@@ -426,6 +417,9 @@
 %override
 param.excludedMetabolites = param.excludedMetabolites | excludedMetabolitesTmp;
 
+if param.printLevel>0
+    disp(param)
+end
 
 %test if the problem is feasible or not
 FBAsolution = optimizeCbModel(model);
@@ -441,17 +435,34 @@
 else
 
     %too time consuming
-    if any(~isfinite(model.S),'all')
-        [I,J]=find(~isfinite(model.S))
-        error('model.S has infinite entries')
-    end
-    if any(~isfinite(model.b))
-        [I,J]=find(~isfinite(model.b))
-        error('model.b has infinite entries')
-    end
-    if any(~isfinite(model.c))
-        [I,J]=find(~isfinite(model.c))
-        error('model.c has infinite entries')
+    if isMATLABReleaseOlderThan('R2022a')
+        if size(model.S,1)*size(model.S,2)<=1e4
+            if any(~isfinite(model.S),'all')
+                [I,J]=find(~isfinite(model.S))
+                error('model.S has non-finite entries')
+            end
+            if any(~isfinite(model.b))
+                [I,J]=find(~isfinite(model.b))
+                error('model.b has non-finite entries')
+            end
+            if any(~isfinite(model.c))
+                [I,J]=find(~isfinite(model.c))
+                error('model.c has non-finite entries')
+            end
+        end
+    else
+        if ~allfinite(model.S)
+            [I,J]=find(~isfinite(model.S))
+            error('model.S has non-finite entries')
+        end
+        if ~allfinite(model.b)
+            [I,J]=find(~isfinite(model.b))
+            error('model.b has non-finite entries')
+        end
+        if ~allfinite(model.c)
+            [I,J]=find(~isfinite(model.c))
+            error('model.c has non-finite entries')
+        end
     end
 
     solution = relaxFBA_cappedL1(model,param);
@@ -495,7 +506,7 @@
                         printConstraints(model,-inf,inf, solution.q>=feasTol,relaxedModel, 0);
                     end
                     if param.relaxedPrintLevel>0 && any(abs(solution.r)>=feasTol)
-                        fprintf('%s\n','The  steady state constraint on this metabolite had to be relaxed:')
+                        fprintf('%s\n','The  steady state constraints on these metabolites had to be relaxed:')
                         disp(model.mets(abs(solution.r)>=feasTol));
                     end
                 end