% 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: 
% 
%
% Inputs: rawRxnEqnString text reaction string as written in FBA odel flat file
%
% Outputs: metabolite specifications for FBA matrix
%
%
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
%
% Notes: 
%
%
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

function [reactants, products, reversible, compartment, failureFlag ] = FBA_EB_rxnEqnParse(rawRxnEqnString)

% initialize outputs
failureFlag = 0;
compartment = '';
reversible = []; % this changes size later, which is not ideal
reactants = {}; % this changes size later, which is not ideal
products = {}; % this changes size later, which is not ideal


% Determine if single compartment rxn
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
idxColon = strfind(rawRxnEqnString, '] : ');
% check for single compartment designation found
if ~isempty(idxColon)
    % compartment found store letter designation
    compartment = rawRxnEqnString( 2 : (idxColon - 1) );
    eqnOnlyRxnString = rawRxnEqnString( (idxColon + 4) : end );
    rxnSingleCompartment = 1;
else
    % no single compartment found, either transport or improperly formatted
    compartment = 't'; %transport
    eqnOnlyRxnString = rawRxnEqnString;
    rxnSingleCompartment = 0;
            
end % if colon separating compartment is found
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 



% Determine reversibility
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

idxReversible = strfind( eqnOnlyRxnString, ' <==> ' );
idxIreversible = strfind( eqnOnlyRxnString, ' ==> ' );

% check that separator was found, store reversibility based on which was found, and
% store reactant/product location for use later
if ~isempty(idxReversible)

    % Reversible reaction separator found
    reversible = 1;
    idxEndReactants = idxReversible;
    idxStartProducts = idxReversible + 5;

elseif ~isempty(idxIreversible)

    % ireversible reaction separator found
    reversible = 0;
    idxEndReactants = idxIreversible;
    idxStartProducts = idxIreversible + 4;

else 

    % no seperator found, failure
    failureFlag = 1;

end % reversible/ireversible separator found
    
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 



% Get reactants
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  
if ~failureFlag
    
    % get string containing only reactants
    reactantsStr = eqnOnlyRxnString( 1 : idxEndReactants );
    
    % reactants are consumed
    signStoich = -1;
    
    % match and sort into coefficients and metabolite name array
    if rxnSingleCompartment
        % single compartment identified earlier
        [ reactants , failureFlag] = metaboliteStoichParse(reactantsStr, signStoich, compartment);
    else
        % no single compartment identified earlier
        [ reactants , failureFlag] = metaboliteStoichParse(reactantsStr, signStoich);
    end % if single compartment identified
        
end % if no failure so far get reactants
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 



% Get products
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
if ~failureFlag
    
    % get string containing only reactants, and find matches for:
    % (number) followed a(space) and then (alphanumeric, underscores, and dashes)
    productsStr = eqnOnlyRxnString( idxStartProducts : end );
    
    % positive stoichiometric coefficents when products are produced
    signStoich = 1;
    
    % match and sort into coefficients and metabolite name array
    if rxnSingleCompartment
        % single compartment identified earlier
        [ products , failureFlag] = metaboliteStoichParse(productsStr, signStoich, compartment);
    else
        % no single compartment identified earlier
        [ products , failureFlag] = metaboliteStoichParse(productsStr, signStoich);
    end % if single compartment identified
    
end % if no failure so far get products
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

end % function FBA_EB_rxnEqnParse



function [ metabolites , failureFlag] = metaboliteStoichParse( oneSideOfRxn, signStoich, singleCompartment )
% gets and parses matches from regular expression (from string input of one side of
% reaction)
% 
% expected format 'stoichCoefficent metaboliteName' or in the case of a transport reaction
% 'stoichCoefficent metaboliteName[compartment]' check of equation parts and sorts them into
% a cell array {'metaboliteName' 'compartment' [stoichCoefficent]} where metabolite name is now
% corrected to include unerscores in place of dashes and is in all caps

% initialize with no failure
failureFlag = 0;

% check for number of inputs
switch nargin
    
    case 3 % single compartment reaction
        % single compartment designated, use for all metabolites
        compartment = singleCompartment;
        if ~isletter(compartment)
            % invalid, fail and add warning about this variable name
            failureFlag = 1;
        end % if compartment invalid
        
        % get metabolite matches for:
        % (number) followed a(space) and then (alphanumeric, underscores, and dashes)
        [~,~,~, metaboliteMatches,~,~, metaboliteLeftovers] = regexp(oneSideOfRxn, '\d+(\.\d+)? (\w|\-)+');

        % number of metabolite matches
        nMetabolites = length(metaboliteMatches);

        if nMetabolites > 0


            % initialize metabolite array
            metabolites{nMetabolites,3} = [];

            % fill looping through metabolite matches
            for iMet = 1 : nMetabolites
                % string for iMet'th metabolite match
                stoichAndMetaboliteString = metaboliteMatches{iMet};
                % find space between stoichiometric coeff and the metabolite name
                idxSpace = strfind( stoichAndMetaboliteString, ' ' );
                % convert string to number
                metaboliteStoich = signStoich * str2double(stoichAndMetaboliteString(1 : idxSpace - 1));
                % replace any dashes in metabolite name with underscore
                metaboliteName = upper(strrep(stoichAndMetaboliteString( (idxSpace + 1) : end), '-', '_'));
                % store
                metabolites(iMet,:) = [{metaboliteName} {compartment} {metaboliteStoich}];

                % failure if metabolite name invalid
                if ~isvarname([ compartment, metaboliteName]) % EWB changed 4/18/2011 to allow metolite names begining with numbers assuming always used as var name preceeded by compartment designation
                    failureFlag = 1;
                end % if metabolite name not valid variable name
                
                % if metabolite has reserved keyword
                if FBA_EB_reservedKeywordCheck(metaboliteName)
                    failureFlag = 1;
                end % if reserved keyword in name

            end % for iMet
            
            % check to make sure nothing potentially left out 
            % (characters other than epected +, space, and seperators in
            if ~isempty( regexp(cat(2,metaboliteLeftovers{:}), '[^+ ]','once') );
                % somehting is left over after regular expresion matching that is not expected
                failureFlag = 1;
            end % if unexpected regular expression leftovers

        else

            % no metabolites found
            metabolites = {};
            failureFlag = 1;

        end % if metabolites found
        
    case 2 % may be transport reaction
        
        % get metabolite matches for:
        % (number) followed a(space) and then (alphanumeric, underscores, and dashes) the
        % ([compartment])
        [~,~,~, metaboliteMatches,~,~, metaboliteLeftovers] = regexp(oneSideOfRxn, '\d+(\.\d+)? (\w|\-)+\[\w\]');

        % number of metabolite matches
        nMetabolites = length(metaboliteMatches);

        if nMetabolites > 0

            % initialize metabolite array
            metabolites{nMetabolites,3} = [];

            % fill looping through metabolite matches
            for iMet = 1 : nMetabolites
                % string for iMet'th metabolite match
                stoichAndMetaboliteString = metaboliteMatches{iMet};
                % find space between stoichiometric coeff and the metabolite name
                idxSpace = strfind( stoichAndMetaboliteString, ' ' );
                % convert string to number
                metaboliteStoich = signStoich * str2double(stoichAndMetaboliteString(1 : idxSpace - 1));
                % replace any dashes in metabolite name with underscore
                metaboliteName = upper(strrep(stoichAndMetaboliteString( (idxSpace + 1) : end - 3), '-', '_'));
                % get compartment & test validity
                compartment = stoichAndMetaboliteString(end - 1);
                if ~isletter(compartment)
                    % invalid, fail and add warning about this variable name
                    failureFlag = 1;
                end % if compartment invalid
                % store
                metabolites(iMet,:) = [{metaboliteName} {compartment} {metaboliteStoich}];

                % failure if metabolite name invalid
                if ~isvarname([ compartment, metaboliteName]) % EWB changed 4/18/2011 to allow metolite names begining with numbers assuming always used as var name preceeded by compartment designation
                    failureFlag = 1;
                end % if metabolite name not valid variable name

                % if metabolite has reserved keyword
                if FBA_EB_reservedKeywordCheck(metaboliteName)
                    failureFlag = 1;
                end % if reserved keyword in name
                
            end % for iMet
            
            % check to make sure nothing potentially left out 
            % (characters other than epected +, space, and seperators in
            if ~isempty( regexp(cat(2,metaboliteLeftovers{:}), '[^+ <>==]','once') );
                % somehting is left over after regular expresion matching that is not expected
                failureFlag = 1;
            end % if unexpected regular expression leftovers

        else

            % no metabolites found
            metabolites = {};
            failureFlag = 1;

        end % if metabolites found
        
    otherwise % too few or too many inputs
        % potetnially throw error...
        
end % switch number of inputs

end % function metaboliteStoichParse