
function [fluxesMaxedProportionally] = FBA_EB_distributeMetabolites(b_vBoundSet, vODEs, S_HV )

% Author: Elsa W. Birch, ewbirch@stanford.edu, Covert Lab - Stanford Bioengineering
% Date: 4/24/2011
% Version: associated with production of the figures for Birch et al. 2012
%
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
%
% Summary: This function calculates the bounding reaction rates for viral
%          synthesis given the metabolite rate vector of host state
%          production and the unlimeted ODE rates using the overall
%          stoichiometric matrix.
% 
%
% Inputs: b_vBoundSet - metabolite rate vector representing host production
%                       of metabolites.
%         vODEs - the unlimited rate 'request' of the viral ODEs
%         S_HV - host viral stoichiometric matrix. This is the combined
%                stoichometric matrix *without* the viral metabolites and
%                *without* any production terms associated with viral reactions.
%                Only consuption terms.
%
% Outputs: fluxesMaxedProportionally - a vector of the maximum rates
%                                      through viral synthesis reactions,
%                                      in the format and order of the
%                                      stoichometric equation.
%
%
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
%
% Notes: 
%
%
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -


% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -


% initialize variable that will eventually be reassigned to output
fluxesWhileRunning = zeros(size(vODEs));

% iteration counter
iCount = 0;

% flux multiplier
ci = 1;
ciMax = 1;

zeroThresh = -1e-5;
logicalExhaustedMetabolites = b_vBoundSet > zeroThresh;
exhaustedStoichRows = S_HV( logicalExhaustedMetabolites, : );
logicalRxnsUsingExhausted = any( (exhaustedStoichRows < 0) , 1);
vODEs( logicalRxnsUsingExhausted ) = zeros(nnz( logicalRxnsUsingExhausted), 1 );

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -


while (ciMax > 0) && (ci > 0) && any(vODEs)

    % Step A
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

    % metabolite request for current fluxes vODEs
    b_vODEs = S_HV * vODEs;
    
    % rounding error fix
    b_vODEs = b_vODEs .* ( b_vODEs < zeroThresh);
    
    % check if all fluxes are within bounds
    allFluxWithinBounds = ~ any( (((b_vODEs) * ciMax) < b_vBoundSet) & ( b_vODEs < zeroThresh ));
    noFluxDrawFromHost = ~ any( (b_vODEs) ); % this is true when all (S_HV * vODEs) = 0
    
    if allFluxWithinBounds || noFluxDrawFromHost
        
        % case when all ( maximum possible as given/remaining of ODEs) fluxes are 
        % currently within metabolite bounds, ci = ciMax and this will be the last
        % itteration
        ci = ciMax;
        
    else
        
        % case when one or more flux is at/above set maximum, so calculate ci
        
        % logical idexing vector with true where ciMax * b_vODEs < b_vBoundSet, these are
        % metabolites for which request by vrius exceeds possible supply
        logicalNonZero = (b_vBoundSet < zeroThresh) & (b_vODEs < zeroThresh);
        logicalODEsRequestExceedsSupply = ( (b_vODEs * ciMax) < b_vBoundSet )  & logicalNonZero;
        if any(logicalODEsRequestExceedsSupply)
            % the largest possible ci for this itteration is the minimum of the fractions of
            % possible/requested metabolites exceeding supply
            cisPossible = b_vBoundSet( logicalODEsRequestExceedsSupply ) ./ ...
                                b_vODEs( logicalODEsRequestExceedsSupply );

            % set ci as 
            ci = min( cisPossible );

            % in the chance that multiple same
            if length(ci) > 1
               ci = ci(1); 
            end
            
        else
            disp('line 91 - look to see why this is getting used?!')
            ci = 0 ;
            
        end
    end % allFluxWithinBounds
    
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
    
    
    % Update output & max
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
    fluxesWhileRunning = fluxesWhileRunning + ( vODEs * ci );
    ciMax = ciMax - ci;

    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
    

    % Step B
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
    
    % remove metabolites alotted to virus durring this iteration from total virus
    % allocation
    b_vBoundSet = b_vBoundSet - (( S_HV * vODEs) * ci );
    
    
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 


    % Step C
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
	
    logicalExhaustedMetabolites = b_vBoundSet > zeroThresh;
    exhaustedStoichRows = S_HV( logicalExhaustedMetabolites, : );
    logicalRxnsUsingExhausted = any( (exhaustedStoichRows < 0) , 1);
    vODEs( logicalRxnsUsingExhausted ) = zeros(nnz( logicalRxnsUsingExhausted), 1 );
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
    
    
    % counter, just for debugging
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
    
    iCount = iCount + 1;
    
    % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

end % while


fluxesMaxedProportionally = fluxesWhileRunning;
