cba_adjust_fluxes

 CBA_ADJUST_FLUXES - Correct a flux mode to make it economical

 v_feasible = cba_adjust_fluxes(v, N, external, cba_constraints, network)

 Adjust given flux distribution v -> economical flux distribution v_feasible

 Input 'network' is only necessary if v contains inactive reactions

0001 function v_feasible = cba_adjust_fluxes(v, N, external, cba_constraints, network)
0002 
0003 % CBA_ADJUST_FLUXES - Correct a flux mode to make it economical
0004 %
0005 % v_feasible = cba_adjust_fluxes(v, N, external, cba_constraints, network)
0006 %
0007 % Adjust given flux distribution v -> economical flux distribution v_feasible
0008 %
0009 % Input 'network' is only necessary if v contains inactive reactions
0010 
0011 if cba_constraints.zv'*v <= 0, 
0012   
0013   warning('Flux distribution v is not beneficial; adjustment is not possible');
0014   v_feasible = []; 
0015 
0016 else
0017 
0018   epsilon    = 10^-8;
0019 
0020   %% replace approximate zeros by zero
0021   [this_v,ind_inactive] = v_exact_zeros(v,N,external,epsilon);
0022   
0023   ind_int = find(external==0); 
0024   N_int   = N(ind_int,:);
0025   [v_act, N_int_act, Es_act, nn_act, cba_constraints_act, ind_act] = cba_reduce_to_active_subnetwork(v, N_int, [], network, cba_constraints);
0026   CC = network_efmtool(nn_act, 'internal', [], cba_constraints_act.zv);
0027   CC(abs(CC)<10^-10) = 0;
0028   C_cba = zeros(length(v),size(CC,2));
0029   C_cba(ind_act,:) = CC;
0030   
0031   % which of the cycles (columns of C_cba) have the same signs as v on their entire support?
0032   ind_violate = find( sign(this_v)' * C_cba == sum(abs(C_cba)) );
0033   C_cba       = C_cba(:,ind_violate);
0034   
0035   while length(C_cba(:)),
0036     fprintf('%d conflicts - ',size(C_cba,2));
0037     clear kappa score
0038 
0039     for it = 1:size(C_cba,2),
0040       gamma                   = C_cba(:,it);
0041       [kappa(it),ind_min(it)] = nanmin(this_v(find(gamma))./gamma(find(gamma)));
0042       score(it)               = kappa(it).^2 * gamma'*gamma;
0043     end
0044     
0045     [score_min,ind]       = min(score);
0046     [this_v,ind_inactive] = v_exact_zeros(this_v - kappa(ind) * C_cba(:,ind),N,external,epsilon);
0047 
0048     %% keep only cycles whose support is inside the support of the updated v
0049     C_cba                 = C_cba(:,sum(abs(C_cba(ind_inactive,:)),1) ==0);
0050 
0051     %% keep only cycles that have the same signs as v on their entire support
0052     ind_violate           = find(sign(this_v)' * C_cba == sum(abs(C_cba)));
0053     C_cba                 = C_cba(:,ind_violate);
0054   end
0055   
0056   v_feasible = this_v;
0057   
0058 end