% 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: Remove futile cycles of metabolites, the potential for which is etermined during parsing
% 
%
% Inputs: 	organism - simulation structure
%			iSim - itteration 
%
% Outputs: 
%
%
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
%
% Notes: 
%
%
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
function [cleanedFluxes] = FBA_EB_cycleRemove(organism, iSim)

cleanedFluxes = organism.simulation.fluxes( :, iSim );

% remove cycles

% loop through cycles
for iCyc = 1 : size( organism.model.potentialCycles , 1 )
    
    idxRxnsCycle = organism.model.potentialCycles{iCyc};

    % get those rxn rows from the stoichiometric matrix
    cycleRxns = organism.model.stoichMatrix(:, idxRxnsCycle );
    
    cycleFluxes = cleanedFluxes( idxRxnsCycle );
    
    cycleNetMet = cycleRxns * cycleFluxes;

    
    nRxnsCycle =  length( idxRxnsCycle );

    % zero out all the cycle fluxes
    cleanedFluxes(idxRxnsCycle) = zeros( nRxnsCycle,1 );
    
    % use firt met to determine net through a single rxn
    idxKeyMet = find( cycleNetMet, 1, 'first' );
    if ~isempty( idxKeyMet );
        
        keyNetMet = cycleNetMet( idxKeyMet);
        
        % cycle through rxns, asigning net flux to them sequentially (eg,
        % if net is less than max of first then only it will have a nonzero
        % flux)
        iRxn = 0;
        while (iRxn < nRxnsCycle) && (abs(keyNetMet) > 0)
            
            % update counter
            iRxn = iRxn + 1;
            
            % stoichiometric coeff of key metabolite in rxn
            stoichRxnKeyMet = organism.model.stoichMatrix(idxKeyMet, idxRxnsCycle(iRxn) );
            
            % flux through that rxn that would result in the net change in
            % key metabolite
            rxnFlux = keyNetMet / stoichRxnKeyMet;
            
            % test for feasibility w/in flux bounds of rxn, and assign
            if rxnFlux > organism.simulation.fluxBoundsUpper( idxRxnsCycle(iRxn), iSim )
                
                % rxn Flux to produce net of key metabolite exceedes rxn
                % bounds, set at bound and assign remaining net key
                % metabolite to subsequent rxn
                
                cleanedFluxes( idxRxnsCycle(iRxn) ) = organism.simulation.fluxBoundsUpper( idxRxnsCycle(iRxn), iSim );
                
                keyNetMet = keyNetMet - stoichRxnKeyMet * cleanedFluxes( idxRxnsCycle(iRxn) );
                
            elseif rxnFlux < organism.simulation.fluxBoundsLower( idxRxnsCycle(iRxn), iSim )
                
                % rxn Flux to produce net of key metabolite exceedes rxn
                % bounds, set at bound and assign remaining net key
                % metabolite to subsequent rxn
                
                cleanedFluxes( idxRxnsCycle(iRxn) ) = organism.simulation.fluxBoundsLower( idxRxnsCycle(iRxn), iSim );
                
                keyNetMet = keyNetMet - stoichRxnKeyMet * cleanedFluxes( idxRxnsCycle(iRxn) );
                
            else
                % rxn Flux to produce net of key metabolite within bounds
                
                cleanedFluxes( idxRxnsCycle(iRxn) ) = rxnFlux;
                
                keyNetMet = keyNetMet - stoichRxnKeyMet * cleanedFluxes( idxRxnsCycle(iRxn) );
                
            end % tests 
            
 
        end % iRxn & keyNetMet 
        
    end % if nonzero net
    
    
end % for iCyc


% % User defined
% % loop through cycles
% for iCyc = 1 : size( organism.model.cycleIdxsUser , 1 )
%     
%     idxRxnsCycle = organism.model.cycleIdxsUser{iCyc};
% %     coeffRxnsCycle = organism.model.cycleCoeffsUser{iCyc};
% 
%     % get those rxn rows from the stoichiometric matrix
%     cycleRxns = organism.model.stoichMatrix(:, idxRxnsCycle );
%     
%     cycleFluxes = cleanedFluxes( idxRxnsCycle );
%     
%     cycleNetMet = cycleRxns * cycleFluxes;
%     find(cycleNetMet)
%     cycleNetMet(find(cycleNetMet))
%     
%     nRxnsCycle =  length( idxRxnsCycle );
% 
%     if all( abs(cycleNetMet) < 1e-7 )        
%         % zero out all the cycle fluxes
%         cleanedFluxes(idxRxnsCycle) = zeros( nRxnsCycle,1 );
%         disp(['cyc remove zeroed fail at',num2str(iSim)])
%     else
%         
%     idxKeyMet = find( cycleNetMet, 1, 'first' );
%     if ~isempty( idxKeyMet );
%         
%         keyNetMet = cycleNetMet( idxKeyMet);
%         
% %         stoichiometric coeff of key metabolite in rxn
%             stoichRxnKeyMet = organism.model.stoichMatrix(idxKeyMet, idxRxnsCycle(1) );
%             
%             % flux through that rxn that would result in the net change in
%             % key metabolite
%             rxnFlux = keyNetMet / stoichRxnKeyMet;
%             
%             % test for feasibility w/in flux bounds of rxn, and assign
%             if ~((rxnFlux < organism.simulation.fluxBoundsLower( idxRxnsCycle(1), iSim )) || (rxnFlux > organism.simulation.fluxBoundsUpper( idxRxnsCycle(1), iSim )))
%                 
%                 cleanedFluxesTentative  = organism.simulation.fluxes( :, iSim );
%                 cleanedFluxesTentative(idxRxnsCycle) = zeros( nRxnsCycle,1 );
%                 cleanedFluxesTentative(idxRxnsCycle(1)) = rxnFlux ;
%                 
%                 cycleNetMetTest = cycleRxns * cleanedFluxesTentative(idxRxnsCycle);
%                 tDif = cycleNetMet - cycleNetMetTest ;
%                 find( tDif )
%                 tDif( find( tDif ) )
%                 testMet = abs( cycleNetMet - cycleNetMetTest ) < 1e-7;
%                 if( all(testMet) )
%                     cleanedFluxes(idxRxnsCycle) = zeros( nRxnsCycle,1 );
%                     cleanedFluxes( idxRxnsCycle(iRxn) ) = rxnFlux;
%                 else
%                     disp(['cyc remove user fail at',num2str(iSim)])
%                 end
%                 
%             end % tests 
%             
%         
%     end % if key met found
%         
%     end
%     
%     
% 
%             
%             
% %     temporary check
%     newRHS = cycleRxns * cleanedFluxes(idxRxnsCycle);
%     oldRHS = cycleRxns * organism.simulation.fluxes( idxRxnsCycle, iSim );
%     checkVal = (newRHS - oldRHS);
%     if ~isempty(find(checkVal))
%         disp(' - - - - -  - - - - -  - - - - -  - - - - -  - - - - -  - - - - -  - - - - - ')
%         c = find(checkVal)
%         checkVal(c)
%         disp(' - - - - -  - - - - -  - - - - -  - - - - -  - - - - -  - - - - -  - - - - - ')
%     end
%  
%     
%     
% end % for iCyc
% % 





