% Author: Elsa W. Birch, ewbirch@stanford.edu, Covert Lab - Stanford Bioengineering
% Date: July 27th 2012
% Version: associated with production of the figures for Birch et al. 2012 
%
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
%
% Summary: Main function to parse in viral model from flat files and functions and set up 
%			the connections required for integrated simulation with FBA
% 
% Inputs: 	viralModelFolder - location of model component fucntions, reaction defintions etc
%			viralModelName - name prefix of viral model
%			organismModelStruct - structure for uninfected FBA organism definition
%			viralParseOptions - structure of parameters and otions that are applied to the aprsing process, 
%								such as the default bounds set for viral reactions
%
% Outputs: 	organismAndVirusModelStruct - infected organism defintion structure for use in integrated simulation
%			viralParseWarnings - cell array for string outputs for warnings that occur durring parsing
%
%
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
%
% Notes: 
%
%
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

function [organismAndVirusModelStruct, viralParseWarnings] = ...
    FBA_EB_viralFileParse( viralModelFolder, viralModelName, organismModelStruct, viralParseOptions )

% in itialize error storage
viralParseWarnings = {};
viralParseErrors = {};
continueParse = 1;

% initialize output
organismAndVirusModelStruct = organismModelStruct;

% store number of host/rxns and mets
organismModelStruct.nHostRxns = organismModelStruct.nRxns;
organismModelStruct.nHostMetabolites = organismModelStruct.nMetabolites;

% file prefix is letters preceeding 'Model' in folder
% [ ~ , viralModelFolderName, ~ ] = fileparts( viralModelFolder );
% viralModelNamePrefix = viralModelFolderName( 1:end-10 ); % 10 because 10 letters in VirusModel
viralModelNamePrefix = viralModelName;
% load in and parse reaction stoichiometry file
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

% first read in just the header row, to determine number of columns to read in
fileIdRxns = fopen( fullfile( viralModelFolder, [viralModelNamePrefix 'VirusStoich.txt'] ));
rawTextReadViralHeaderFind = textscan(fileIdRxns, '%s', 'delimiter', '\n');
fclose(fileIdRxns); % close file
rawTextReadViralHeader = rawTextReadViralHeaderFind{1};
rawTextReadViralHeader = rawTextReadViralHeader{1};

% 
rxnIDStrIdx = strfind(rawTextReadViralHeader,'rxnID');
inputToHostFlagStrIdx = strfind(rawTextReadViralHeader,'inputToHost');
rxnSepStrIdx = strfind(rawTextReadViralHeader,'<==>');
% check that column identifiers are found in correct order, and devide reactant/product raw strings
if isempty(rxnIDStrIdx) || isempty(inputToHostFlagStrIdx) || isempty(rxnSepStrIdx) 
    
    % one of eqn/ODE identifiers not found
    viralParseErrors = [viralParseErrors; ...
                       {'column identifier ''rxnID'', ''inputToHost'', or ''<==>'' not found '}];
    continueParse = 0;
                   
elseif ~ (( rxnIDStrIdx < inputToHostFlagStrIdx ) && ( inputToHostFlagStrIdx < rxnSepStrIdx ))
    
    % identifiers in wrong order
    viralParseErrors = [viralParseErrors; ...
                       {'column identifier ''rxnID'', ''inputToHost'', and ''<==>'' occur in wrong order '}];
    continueParse = 0;
    
else
    % correct
    rawTextReadViralHeaderReactants = rawTextReadViralHeader( inputToHostFlagStrIdx + 11 : rxnSepStrIdx - 1 );
    rawTextReadViralHeaderProducts = rawTextReadViralHeader( rxnSepStrIdx + 4 : end );
    
    % convert dashes to underscores
    rawTextReadViralHeaderReactants = strrep( rawTextReadViralHeaderReactants, '-', '_');
    rawTextReadViralHeaderProducts = strrep( rawTextReadViralHeaderProducts, '-', '_');
    
end % if, check for column identifiers

% check number of reactants/products
if continueParse
    
    % get metabolite designations
    [~, ~, ~, productMatchStr, ~, ~, productSplitStr] = ...
        regexp(rawTextReadViralHeaderProducts,'\w+\[[a-zA-Z]\]');
    [~, ~, ~, reactantMatchStr, ~, ~, reactantSplitStr] = ...
        regexp(rawTextReadViralHeaderReactants,'\w+\[[a-zA-Z]\]');
    
    % test for unexpected charecters
    [~, ~, ~, unexpctedCharacters, ~, ~, ~] = ...
        regexp(cat(2,productSplitStr{:},reactantSplitStr{:}),'\S');
    if ~isempty( unexpctedCharacters )
        
        % unexpected characters found
        viralParseErrors = [viralParseErrors; ...
                       {['unexpected characters found in reactant or header columns ',[unexpctedCharacters{:}]]}];
        continueParse = 0;
        
        % throw errors
        viralParseErrors = [viralParseErrors];
        if ~isempty(viralParseErrors)
            for iEr = 1 : length(viralParseErrors)
                disp(sprintf('%s', viralParseErrors{iEr} ))
            end % iEr
            error('viral parsing failed, see above')
        end % if errors to throw
        
    else
        
        % get number of each products and reactatants
        nReactantMatches = length(reactantMatchStr);
        nProductMatches = length(productMatchStr);
        
    end % if ~isempty( unexpctedCharacters )
    
    
end %if, check number of reactants/products

% parse rest of viral stoich file
if continueParse % parse rest of viral stoich file
    
    % create expected format string for text import
    reactantFormatString = repmat('%f ',1 ,nReactantMatches);
    productFormatString = repmat('%f ',1 ,nProductMatches);
    parseFormat = ['%s %d ',reactantFormatString,'%s ',productFormatString];
    
    % import text
    fileIdRxns = fopen( fullfile( viralModelFolder,[viralModelNamePrefix 'VirusStoich.txt'] ));
    headerRowsRxns = 1; % # ignored header rows
    rawTextReadViralStoich = textscan(fileIdRxns, parseFormat, 'Delimiter', '\t','HeaderLines', headerRowsRxns );
    fclose(fileIdRxns); % close file
    
end % if continueParse % parse rest of viral stoich file 

% parse viral rxn metabolites
if continueParse
    
    % initializations
    viralReactantIdxTranslation = zeros(nReactantMatches,1);
    viralProductIdxTranslation = zeros(nProductMatches,1);
    
    % add viral as valid compartment
    organismModelStruct.compartmentList = [organismModelStruct.compartmentList , {'v'} ];
    
    % add metabolites
    for iReact = 1 : nReactantMatches
        
        % metabolite name
        metaboliteMatch = reactantMatchStr{iReact};
        [metaboliteIdxInModel, organismModelStruct, viralReactErrors] = ...
            metaboliteParse(metaboliteMatch,organismModelStruct);
        viralReactantIdxTranslation( iReact ) = metaboliteIdxInModel;
            
        viralParseErrors = [viralParseErrors; viralReactErrors];
    end % for iReact
    
    for iProd = 1 : nProductMatches
        
        % metabolite name
        metaboliteMatch = productMatchStr{iProd};
        [metaboliteIdxInModel, organismModelStruct, viralProdErrors] = ...
            metaboliteParse(metaboliteMatch,organismModelStruct);
        viralProductIdxTranslation( iProd ) = metaboliteIdxInModel;
        
        viralParseErrors = [viralParseErrors; viralProdErrors];
            
    end % for iProd
    
    % trow errors collected
%     viralParseErrors = [viralParseErrors; viralProdErrors; viralReactErrors];
    if ~isempty(viralParseErrors)
        for iEr = 1 : length(viralParseErrors)
            disp(sprintf('%s', viralParseErrors{iEr} ))
        end % iEr
        error('viral parsing failed, see above')
    end % if errors to throw


end % if continueParse, viral rxn metabolites

% parse viral rxns
if continueParse
    
    % get reactions from text intake
    rawTextReadViralStoichRxns = rawTextReadViralStoich{ 1 };
    rawTextReadViralStoichInToHost = rawTextReadViralStoich{ 2 };
    rawTextReadViralStoichReversibility = rawTextReadViralStoich{ 3 + nReactantMatches };
    nViralRxns = length( rawTextReadViralStoichRxns );
    
    % initialize idx translation (for ODES)
    organismModelStruct.viralRxnODETranslation = zeros( nViralRxns, 2);
    
    % store number of viral rxns and mets
    organismModelStruct.nViralRxns = nViralRxns;
    organismModelStruct.nViralMetabolites = organismModelStruct.nMetabolites ...
                                                - organismModelStruct.nHostMetabolites;
    
    % loop through rxns
    for iViralRxn = 1 : nViralRxns
        
        viralRxnName = upper( strrep( rawTextReadViralStoichRxns{ iViralRxn }, '-', '_' ) );
        viralRxnNameWithCompartment = ['v', viralRxnName ];
        
        % test for rxn name validity
        if isvarname(viralRxnNameWithCompartment)
            
            % begin to update model
            organismModelStruct.nRxns = organismModelStruct.nRxns + 1;
            organismModelStruct.rxnList{organismModelStruct.nRxns} = viralRxnNameWithCompartment ;
            organismModelStruct.rxnCompartmentList{organismModelStruct.nRxns} = 'v';
            organismModelStruct.rxns.(viralRxnNameWithCompartment).idxRxn = organismModelStruct.nRxns;
            organismModelStruct.rxns.(viralRxnNameWithCompartment).inToHost = (rawTextReadViralStoichInToHost( iViralRxn ) == 1);
            organismModelStruct.inToHost( organismModelStruct.nRxns ) = organismModelStruct.rxns.(viralRxnNameWithCompartment).inToHost;
            organismModelStruct.rxns.(viralRxnNameWithCompartment).metabolicClass = 'viralReactions';
            organismModelStruct.rxns.(viralRxnNameWithCompartment).metabolicSubClass = 'viralRactions';
            
            % determine and store reversibility
            if ~isempty( strfind( rawTextReadViralStoichReversibility{iViralRxn}, '<==>' ) )
                
                % rxn is reversible
                organismModelStruct.reversibility(organismModelStruct.nRxns) = 1;
                organismModelStruct.rxns.(viralRxnNameWithCompartment).reversible = 1;
                
                % set default flux bounds
                organismModelStruct.fluxBoundsDefault( organismModelStruct.nRxns , : ) = viralParseOptions.fluxBoundsDefault;
                
            elseif ~isempty( strfind( rawTextReadViralStoichReversibility{iViralRxn}, '==>' ) )
                
                % rxn is irreversible
                organismModelStruct.reversibility(organismModelStruct.nRxns) = 0;
                organismModelStruct.rxns.(viralRxnNameWithCompartment).reversible = 0;
                
                % set default flux bounds
                organismModelStruct.fluxBoundsDefault( organismModelStruct.nRxns , : ) = [0 viralParseOptions.fluxBoundsDefault(2)];
                
            else
                
                % proper reversibility string not found
                viralParseErrors = [viralParseErrors; ...
                       {['reversibility string, ',rawTextReadViralStoichReversibility{iViralRxn},', is not either ''<==>'' or ''==>'' ']}];
                
            end % if, reversibility check
            
            % update reactant stoich
            for iReact = 1 : nReactantMatches

                stoichCoeff = -1 * rawTextReadViralStoich{ 2 + iReact }(iViralRxn);
                organismModelStruct.stoichMatrix( viralReactantIdxTranslation(iReact) ,...
                    organismModelStruct.nRxns ) = stoichCoeff;
            end % for iReact
            
            % update product stoich
            for iProd = 1 : nProductMatches

                stoichCoeff = rawTextReadViralStoich{ 3 + nReactantMatches + iProd }(iViralRxn);
                organismModelStruct.stoichMatrix( viralProductIdxTranslation(iProd) ,...
                    organismModelStruct.nRxns ) = organismModelStruct.stoichMatrix( viralProductIdxTranslation(iProd) ,...
                    organismModelStruct.nRxns ) + stoichCoeff; % + is in order to not overwrite reactant stoich

            end % for iProd
            
            % set rule to true, will actually be regulated by ODE
            organismModelStruct.rxnRegulationRulesEvalStrings{organismModelStruct.nRxns} = 'true';
            
        else
            % not a valid variable name
            viralParseErrors = [viralParseErrors; ...
                       {['viral rxn name, ',viralRxnNameWithCompartment,', is not valid']}];
            continueParse = 0;
        end % if isvarname(viralRxnNameWithCompartment)
        
    end % for iViralRxn
    
    % throw errors
    if ~isempty(viralParseErrors)
        for iEr = 1 : length(viralParseErrors)
            disp(sprintf('%s', viralParseErrors{iEr} ))
        end % iEr
        error('viral parsing failed, see above')
    end % if errors to throw
    
end % % parse viral rxns

% parse file holding viral ODE parameters and check for ODE function 
if continueParse
    
    % check for the existance of viral ODE handle
    viralODEname = [viralModelNamePrefix,'VirusODEs'];
    % add to path to test
    
    if exist( [viralModelNamePrefix 'VirusODEs'] , 'file') == 2
        
        % if needed call to setup additional viral arguments
        viralAdditionalArgumentFunctionName =  [viralModelNamePrefix 'AdditionalVarargSetup.m'] ;
        if ( exist(viralAdditionalArgumentFunctionName, 'file') == 2 )
            viralAdditionalArgumentHandle = eval( ['@',viralModelNamePrefix 'AdditionalVarargSetup'] );
            organismModelStruct.viralVarargin = viralAdditionalArgumentHandle( [] );
        else
            organismModelStruct.viralVarargin = {};
        end
        
        % determine if each step varargin function exists
        viralAdditionalArgumentEachRunFunctionName =  [viralModelNamePrefix 'VarargEachRun.m'] ;
        if ( exist(viralAdditionalArgumentEachRunFunctionName, 'file') == 2 )
            organismModelStruct.viralEachRunVararginHandle = eval( ['@',viralModelNamePrefix 'VarargEachRun'] );
            organismModelStruct.viralEachRunFlag = true;
        else
            organismModelStruct.viralEachRunVararginHandle = {};
            organismModelStruct.viralEachRunFlag = false;
        end
        
         % determine if initializer function exists
        viralInitializerFunctionName =  [viralModelNamePrefix, 'VirusInitializer.m'] ;
        if ( exist(viralInitializerFunctionName, 'file') == 2 )
            organismModelStruct.viralInitializerHandle = eval( ['@',viralModelNamePrefix, 'VirusInitializer'] );
            organismModelStruct.viralInitializerFlag = true;
        else
            organismModelStruct.viralInitializerHandle = {};
            organismModelStruct.viralInitializerFlag = false;
        end
        
        checkViralODEs = 1;
        
        % try evaluating handle to check for errors
        try
            
            % save handle to viral ODEs
            handleString = ['@',viralODEname];
            handleViralODEs = eval( handleString );
            organismModelStruct.handleViralODEs = handleViralODEs;
            
            % set up viral ODE translation
            % . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 
            translationIDs =  handleViralODEs([],[],'identifiers',organismModelStruct.viralVarargin,{} );
%             organismModelStruct.translationIDs = translationIDs ;
            organismModelStruct.nViralODEs = size( translationIDs, 1 );
            
            %initialize translation counters
            nRxnTrans = 0;
            nMetTrans = 0;
            
            % set ODE translation
            for iTrans = 1 : organismModelStruct.nViralODEs
                
                % for each ODE, determine if metabolite correspondance is being specified
                if ~isempty(translationIDs{ iTrans, 1}) 
                    
                    metaboliteID = ['v', upper( strrep(translationIDs{ iTrans, 1}, '-','_') ) ];
                    
                    % if met is found in fba model, save tranlsation
                    if isfield( organismModelStruct.metabolites ,metaboliteID)
                        
                        nMetTrans = nMetTrans + 1;
                        fbaMetIdx = organismModelStruct.metabolites.(metaboliteID).idxMetabolite;
                        organismModelStruct.viralMetaboliteODETranslation( nMetTrans , : ) = [fbaMetIdx, iTrans];
                        
                    else
                        viralParseErrors = [viralParseErrors; ...
                            {['viral ODEs metaboliteID, ', metaboliteID ,', is not found in FBA model']}];
                    end % if isfield
                end
                
                % for each ODE, determine if rxn correspondance is being specified
                if ~isempty(translationIDs{ iTrans, 2}) 
                    
                    rxnID = ['v', upper( strrep(translationIDs{ iTrans, 2}, '-','_') ) ];
                    
                    % if rxn is found in fba model, save tranlsation
                    if isfield( organismModelStruct.rxns ,rxnID)
                        nRxnTrans = nRxnTrans + 1;
                        fbaRxnIdx = organismModelStruct.rxns.(rxnID).idxRxn;
                        organismModelStruct.viralRxnODETranslation( nRxnTrans , : ) = [fbaRxnIdx, iTrans];
                        
                    else
                        viralParseErrors = [viralParseErrors; ...
                            {['viral ODEs rxnID, ',rxnID,', is not found in FBA model']}];
                    end % if isfield
                end
                
                
                if isempty(translationIDs{ iTrans, 2}) && isempty(translationIDs{ iTrans, 1})
                    
                    % case when ODE is associated with neither a rxn or matabolite
                    % note, for now nothing is done here, may add warnign/error
                    
                end % if no rxn or metabolite specified
                
                
            end % for iTrans
            
            logicalIdxODEsWithTranslation = organismModelStruct.viralRxnODETranslation(:,1) > 0;
            organismModelStruct.viralRxnODETranslation = organismModelStruct.viralRxnODETranslation(logicalIdxODEsWithTranslation,:);
            % . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 
            
            
            % try regular evaluation of ODEs
            % . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 
            viralODELimitsValues = zeros(organismModelStruct.nViralODEs, 1);
            viralODElimitsFlags = false(organismModelStruct.nViralODEs, 1);
            viralODElimits = {viralODElimitsFlags,viralODELimitsValues};
            junk = handleViralODEs(0, zeros( organismModelStruct.nViralODEs, 1 ) , organismModelStruct.viralVarargin, viralODElimits );
        	% . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 
            
        catch viralODEError
            
            % viral ODE function error
            viralParseErrors = [viralParseErrors; ...
                       {['attempt to evaluate ',viralODEname,'( t, C, param, limits, varargin ) produced error']},];
            continueParse = 0;
            
        end % try
        
    else
        
        % viral ODE function not found
        viralParseErrors = [viralParseErrors; ...
                       {['viral ODEs, ',viralODEname,', is not found']}];
        continueParse = 0;
        
    end
    
    
        % throw errors
    if ~isempty(viralParseErrors)
        for iEr = 1 : length(viralParseErrors)
            viralParseErrors{iEr}
            disp(sprintf('%s', viralParseErrors{iEr} ))
        end % iEr
        
        if exist('viralODEError','var')
           rethrow(viralODEError) 
        end
        
        error('viral parsing failed, see above')
    end % if errors to throw
    
    
end % if continueParse viral parameters

% save if no errors so far
if continueParse
    
    
    % find viral-host metabolite connections
    S_HV = zeros( organismModelStruct.nMetabolites, organismModelStruct.nRxns );
    S_HV(  1 : organismModelStruct.nHostMetabolites , (organismModelStruct.nHostRxns + 1) : end ) = ...
        organismModelStruct.stoichMatrix(  1 : organismModelStruct.nHostMetabolites , (organismModelStruct.nHostRxns + 1) : end );
    organismModelStruct.hostViralMetaboliteConnectionsIdx = find( any( S_HV , 2) );
    
    
    % set output now that structure has been modified with no errors
    organismAndVirusModelStruct = organismModelStruct;
    
end % continueParse, save viral model for output

organismAndVirusModelStruct = orderfields(organismAndVirusModelStruct);

end % main function


% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
function [metaboliteIdxInModel, organismModelStruct, viralParseErrors] = metaboliteParse(metaboliteMatch,organismModelStruct)

%initialize error handeling
viralParseErrors = {};
continueParse = 1;

% initialze output
metaboliteIdxInModel = 0;

% separate name and compartment
metaboliteName = upper( strrep( metaboliteMatch( 1: end - 3 ), '-', '_' ) );
metaboliteCompartment = lower(metaboliteMatch( end - 1 ));

% check that compartment is valid
if ~any( strcmp(metaboliteCompartment,organismModelStruct.compartmentList) )
    viralParseErrors = [viralParseErrors; ...
               {['unexpected compartments designation: ',metaboliteCompartment]}];
    continueParse = 0;
end % if invalid compartment

%
metaboliteNameWithCompartment = [metaboliteCompartment,metaboliteName];

% check if already exists in organism
if isfield(organismModelStruct.metabolites, metaboliteNameWithCompartment) && continueParse

    % found in model, set translation info
    metaboliteIdxInModel = organismModelStruct.metabolites.(metaboliteNameWithCompartment).idxMetabolite;
    
else
    
    % not found in model, check if virus
    if strcmp('v',metaboliteCompartment) && continueParse

        % check that compartment+name is valid matlab variable
        if isvarname(metaboliteNameWithCompartment)
            
            % add to model, metabolites
            organismModelStruct.nMetabolites = organismModelStruct.nMetabolites + 1;
            organismModelStruct.metaboliteList{ organismModelStruct.nMetabolites } ...
                                                        = metaboliteNameWithCompartment;
            organismModelStruct.metaboliteCompartmentList{ organismModelStruct.nMetabolites } ...
                                                        = metaboliteCompartment;                                      
            organismModelStruct.metabolites.(metaboliteNameWithCompartment).idxMetabolite ...
                                                        = organismModelStruct.nMetabolites;
            organismModelStruct.metabolites.(metaboliteNameWithCompartment).rxnsUsingMet ...
                                                        = {};
            % add to model stoich matrix, only zeros at this time
            organismModelStruct.stoichMatrix(organismModelStruct.nMetabolites,organismModelStruct.nRxns) = 0;                                     
            
            % set translation info
            metaboliteIdxInModel = organismModelStruct.metabolites.(metaboliteNameWithCompartment).idxMetabolite;
            
        else 
            
            % cannot add because invalid variable name
            viralParseErrors = [viralParseErrors; ...
               {['metabolite, ',metaboliteNameWithCompartment, ', is designated viral, but cannot be added because name is not valid matlab variable: ']}];

            
        end % if valid variable name
        
        
    else
        
        %found niether in model nor designated virus
        viralParseErrors = [viralParseErrors; ...
               {['metabolite is niether viral nor found in model: ',metaboliteNameWithCompartment]}];
    
    end % is v (virus) metabolite

end % if already in model


end % function metaboliteParse
      
