cba_check

 CBA_CHECK - Check flux mode, chem. pot. differences, and spec. flux costs for feasibility

 isFeasible = cba_check(N, external, cba_constraints, cba_options, v, dmu, hv)

 Check whether flux mode v, chemical potentials dmu, and specific flux costs hv
 satisfy the thermodynamic and economic constraints

0001 function isFeasible = cba_check(N, external, cba_constraints, cba_options, v, dmu, hv)
0002 
0003 % CBA_CHECK - Check flux mode, chem. pot. differences, and spec. flux costs for feasibility
0004 %
0005 % isFeasible = cba_check(N, external, cba_constraints, cba_options, v, dmu, hv)
0006 %
0007 % Check whether flux mode v, chemical potentials dmu, and specific flux costs hv
0008 % satisfy the thermodynamic and economic constraints
0009 
0010 % !!!
0011 % care about zero fluxes ???
0012 
0013 %  ------------------------------------------------------------
0014 % initialisation
0015 
0016 %rand('state',cba_options.seed);
0017 
0018 [nm,nr]      = size(N);
0019 
0020 if ~isfield(cba_constraints,'dmu_min'), cba_constraints.dmu_min = -10 *ones(nr,1); end
0021 if ~isfield(cba_constraints,'dmu_max'), cba_constraints.dmu_max =  10 *ones(nr,1);; end
0022 if ~isfield(cba_constraints,'y_min'),   cba_constraints.y_min   = -10 *ones(nr,1);; end
0023 if ~isfield(cba_constraints,'y_max'),   cba_constraints.y_max   =  10 *ones(nr,1);; end
0024 
0025 cba_constraints.v_min(find(cba_constraints.v_sign== 1)) = 0;
0026 cba_constraints.v_max(find(cba_constraints.v_sign==-1)) = 0;
0027 
0028 zv           = cba_constraints.zv;
0029 ind_v_fix    = find(isfinite(cba_constraints.v_fix));
0030 ind_extsigns = find(isfinite(cba_constraints.ext_sign));
0031 n_int        = sum(external==0); 
0032 N_int        = N(find(external==0),:);
0033 Ktot         = analyse_N(N);
0034 Kint         = analyse_N(N(find(external==0),:));
0035 nKint        = size(Kint,2);
0036 nKtot        = size(Ktot,2);
0037 my_eye       = eye(nr);
0038 
0039 weights      = ones(3*nr,1);     % weights in regularisation term
0040 
0041 
0042 % ------------------------------------------------------------
0043 % vector x contains all v, dmu, and hv
0044 % indices within vector:
0045 
0046 ind_v   =          1:nr;
0047 ind_dmu =    nr + (1:nr);
0048 ind_hv  = 2* nr + (1:nr);
0049 
0050 % ------------------------------------------------------
0051 % cba_constraints:
0052 %
0053 % limits           : cba_constraints.v_min   <=  x(ind_v)   <= cba_constraints.v_max
0054 %                  : cba_constraints.dmu_min <=  x(ind_dmu) <= cba_constraints.dmu_max
0055 %                  : cba_constraints.y_min   <=  x(ind_hv)  <= cba_constraints.y_max
0056 %
0057 % fixed values     : x(ind_v_fix)         = cba_constraints.v_fix(ind_v_fix)
0058 % stationary fluxes: N_int * x(ind_v )    = 0
0059 % dmu constraint   : Ktot' * x(ind_hv)    = 0
0060 % hv constraint    : Kint' * x(ind_hv)    = Kint' * zv
0061 %
0062 % feasible signs   : diag(x(ind_v)) * N' * x(ind_dmu) <= - epsilon < 0   where v == 0
0063 %                  : diag(x(ind_v)) * N' * x(ind_hv)  >=   epsilon > 0   where v == 0
0064 %
0065 % regularisation   : min != sum(x.^2);
0066 
0067 
0068 % ------------------------------------------------------
0069 % rewrite this in the form
0070 %
0071 % Aeq * x  = beq
0072 %   A * x <= b
0073 %       C  = cbaf2(x,N)   (for signs)
0074 %     min != cbaf1(x)     (regularisation)
0075 
0076 
0077 xmin = [cba_constraints.v_min; cba_constraints.dmu_min; cba_constraints.y_min];
0078 xmax = [cba_constraints.v_max; cba_constraints.dmu_max; cba_constraints.y_max];
0079 
0080 Aeq  = double([my_eye(ind_v_fix,:),  zeros(length(ind_v_fix),2*nr)      ; ...
0081                N_int,                 zeros(n_int,2*nr)                ; ...
0082                zeros(nKtot,nr),       Ktot',          zeros(nKtot,nr)  ; ...
0083                zeros(nKint,2*nr),                     Kint'             ]);
0084 
0085 beq  = double([cba_constraints.v_fix(ind_v_fix);... 
0086                zeros(n_int,1)           ;...
0087                zeros(nKtot,1)    ;...
0088                Kint' * zv ]);
0089 
0090 A = double([-diag(cba_constraints.ext_sign(ind_extsigns)) * N(ind_extsigns,:), ...
0091             zeros(length(ind_extsigns),2*nr)]);
0092 
0093 b = zeros(length(ind_extsigns),1);
0094 
0095 
0096 % ------------------------------------------------------------------------
0097 
0098 x = [v; dmu; hv];
0099 
0100 isFeasible = 1;
0101 
0102 ppp.N = full(N);
0103 ppp.ind_v   = ind_v  ;
0104 ppp.ind_dmu = ind_dmu;
0105 ppp.ind_hv  = ind_hv ;
0106 ppp.epsilon = 0;
0107 
0108 if sum(x < xmin) + sum(x > xmax),          warning('Limit constraints violated');      isFeasible = 0; end 
0109 if norm( Aeq * x - beq)/length(beq)>10^-5, warning('Equality constraints violated');   isFeasible = 0; end 
0110 if sum(A * x > b),                         warning('Inequality constraints violated'); isFeasible = 0; end 
0111 if sum(cbaf2(x,ppp) > 0),                  warning('Sign constraints violated');       isFeasible = 0; end