function out = bns_pe(varargin)

% used to parse eventlog files for 1 or more runs
% all should have the same plot_info (same number of plots, same compounds in
% each plot) and the same run_info (same time points)

% If 2 arguments, 1st must be full path to output directory and 2nd must be
% the setsname (currently called output_filename)

% If 0 arguments, it will get output directory and output_filename from default
% directories.

% Modified to work with new output format   v3.0
%                           Slava  02/07

% Modified to work with agregated number of reaction events  v3.1 
%                           Slava  08/07


if nargin == 2
    output_dir = varargin{1};
    data_partname = [varargin{2} '_data'];
    evlog_partname = [varargin{2} '_eventlog'];
    rundata_file = [varargin{2} '_rundata.mat'];
    exec_file = [varargin{2} '_timesexecuted'];
    output_filename = varargin{2};
elseif nargin == 0
    default_directories
    
    if isempty(output_directory)
        output_directory = 'output';
        output_dir = fullfile(models_dir,specific_model_dir,parameters_dir,output_directory);
    else
        output_dir = output_directory;
    end
    data_partname = [output_filename '_data'];
    evlog_partname = [output_filename '_eventlog'];
    rundata_file = [output_filename '_rundata.mat'];
    exec_file = [varargin{2} '_timesexecuted'];
else
    disp('Incorrect number of arguments.  Zero arguments means to use the output file specified in default_directories.m.')    
    disp('If two arguments, they should (1) be the full pathname of the output file and (2) the output_filename.')    
    return
end

% Delete old parsed data with current output_filename. Otherwise there is confusion about which files
% were from which runs.  Use a new output_filename if you don't want to lose data.
parsedname = [output_filename '_parsed'];
delete(fullfile(output_dir, [parsedname '*.mat']));

datafile_list = dir(fullfile(output_dir, [data_partname '*']));
num_datafiles = size(datafile_list,1);

global stopwait
stopwait=false;
%hwait = waitbar(0.0,'Parsing Event Log Files...','CreateCancelBtn','global stopwait; stopwait = true;');
%set(hwait,'Name','Parsing Event Log Progress');


if exist(fullfile(output_dir, rundata_file),'file')
    load(fullfile(output_dir, rundata_file), 'plot_info', 'run_info', 'rxn_grp', 'cmpd_info')
else
    error([fullfile(output_dir, rundata_file) ' does not exist'])
end

monitor_cmpds_gid = plot_info.monitor_cmpds_gid;
eventlog_maxelements = plot_info.eventlog_maxelements;
time_start = run_info.time_start;
num_cmpds_monitored = length(monitor_cmpds_gid);

numgrps = length(rxn_grp);
run_id = 1;

stoich_data = [];
cmpds_affected = [];
monlist_idx = [];
num_cmpds_affected = [];
for ngrp = 1:numgrps
    rxgp = rxn_grp{ngrp};
    numrxns(ngrp) = length(rxgp.rxn);
    max_gid = max(rxgp.global_comp_gid);
    for i = 1:numrxns(ngrp)
        num_cmpds_affected(i,ngrp) = 0;
        rxn = rxgp.rxn{i};
        num_subs = size(rxn.substrates,1);
        for j = 1:num_subs
            id = rxn.substrates(j,1);
            this_gid = rxgp.global_comp_gid(id);
            mon_this = find(this_gid == monitor_cmpds_gid);
            if ~isempty(mon_this)
                num_cmpds_affected(i,ngrp) = num_cmpds_affected(i,ngrp) + 1;
                monlist_idx(i,num_cmpds_affected(i,ngrp),ngrp) = mon_this;
                cmpds_affected(i,num_cmpds_affected(i,ngrp),ngrp) = this_gid;
                switch rxn.type
                    case {'binding','unbinding','uni-uni','uni-multi','homo-ltor','unique-bi-bi'}
                        stoich_data(i,num_cmpds_affected(i,ngrp),ngrp) = -1;
                    case {'n-order','transv1','genmm','undef'}
                        stoich_data(i,num_cmpds_affected(i,ngrp),ngrp) = -rxn.substrates(j,2);
                end
            end
        end
        if isfield(rxn, 'products')
            num_prods = size(rxn.products,1);
            for j = 1:num_prods
                id = rxn.products(j,1);
                this_gid = rxgp.global_comp_gid(id);
                mon_this = find(this_gid == monitor_cmpds_gid);
                if ~isempty(mon_this)
                    cmpd_ex = 0;   
                    for i_ex=1:num_cmpds_affected(i,ngrp)
                       if this_gid == cmpds_affected(i,i_ex,ngrp) 
                         cmpd_ex = i_ex;
                         break
                       end
                    end
                  if (cmpd_ex > 0)
                      switch rxn.type
                        case {'binding','unbinding','uni-uni','uni-multi','homo-ltor','unique-bi-bi'}
                            stoich_data(i,cmpd_ex,ngrp) = stoich_data(i,cmpd_ex,ngrp) + 1;
                        case {'n-order','transv1','genmm','undef'}
                            stoich_data(i,cmpd_ex,ngrp) = stoich_data(i,cmpd_ex,ngrp) + rxn.products(j,2);
                      end
                  else
                    num_cmpds_affected(i,ngrp) = num_cmpds_affected(i,ngrp) + 1;
                    monlist_idx(i,num_cmpds_affected(i,ngrp),ngrp) = mon_this;
                    cmpds_affected(i,num_cmpds_affected(i,ngrp),ngrp) = this_gid;
                    switch rxn.type
                        case {'binding','unbinding','uni-uni','uni-multi','homo-ltor','unique-bi-bi'}
                            stoich_data(i,num_cmpds_affected(i,ngrp),ngrp) = 1;
                        case {'n-order','transv1','genmm','undef'}
                            stoich_data(i,num_cmpds_affected(i,ngrp),ngrp) = rxn.products(j,2);
                    end
                  end
                end
            end
        end
        

    end
end

execfile_list = dir(fullfile(output_dir, [exec_file '*']));
num_execfiles = size(execfile_list,1);        
times_rxns_executed = [];

for nf = 1:num_execfiles
    if exist(fullfile(output_dir, execfile_list(nf).name),'file')
        times_rxns_executed = read_timesexecuted(fullfile(output_dir, execfile_list(nf).name), times_rxns_executed);
    else
        error([fullfile(output_dir, execfile_list(nf).name) ' does not exist'])
    end
end

evlog_part2name = [evlog_partname '_' num2str(run_id)];
ev_log_name = [evlog_part2name '_' num2str(1)];
if exist(fullfile(output_dir, ev_log_name),'file')
    event_log = read_logdata(fullfile(output_dir, ev_log_name));

if size(event_log,2) == 3
    grpnum_exists = 1;
else
    grpnum_exists = 0;
end

for run_id = 1:num_datafiles
 
    if stopwait 
        delete(hwait)
        clear hwait
        break; 
    end
%    waitbar(nrun/num_datafiles,hwait,['Run Number ' num2str(nrun)]);
    
    data_time = [];
    data_cmpds = [];
    temp_data_time = [];
    temp_data_cmpds = [];
    for i = 1:num_cmpds_monitored
        data_time{i} = [];
        data_cmpds{i} = [];
        each_cmpd_events(i) = 0;
    end

    cmpd_vals = cmpd_info.orig_cmpd_vals;
   
    evlog_part2name = [evlog_partname '_' num2str(run_id)];
    elfile_list = dir(fullfile(output_dir, [evlog_part2name '_*']));
    num_elfiles = size(elfile_list,1);
    
    temp_rxntimes = [];
    max_rxnevents = 0;
    max_rxnum = 0;
    tot_events = 0;
    for ngrp = 1:numgrps
        for nrx = 1:numrxns(ngrp)
            rxn_eventtimes{nrx,ngrp} = [];
            n_event_rxn(nrx,ngrp) = 0;
            tot_events = tot_events + times_rxns_executed{ngrp,run_id}(nrx);
            for ncmpd = 1:num_cmpds_affected(nrx,ngrp)
                c_idx = monlist_idx(nrx,ncmpd,ngrp);
                each_cmpd_events(c_idx) = each_cmpd_events(c_idx) + times_rxns_executed{ngrp, run_id}(nrx);
            end
        end
        max_rxnevents = max([max_rxnevents max(times_rxns_executed{ngrp, run_id})]);
        max_rxnum = max([max_rxnum numrxns(ngrp)]);
    end
    tot_cmpd_events = sum(each_cmpd_events);
    max_cmpd_events = max(each_cmpd_events);
    
    if tot_events > eventlog_maxelements
        add_elsize = eventlog_maxelements;
    else
        add_elsize = tot_events;
    end
    if grpnum_exists
        add_elsize = add_elsize*3;
    else
        add_elsize = add_elsize*2;
    end
    
    tot_elements_needed = tot_events + 2*tot_cmpd_events + add_elsize;   % for rxn_eventtimes, data_time, data_cmpds, and eventlog.
    
    % check if internal matrices are small enough for fast parsing
    if max_rxnum*max_rxnevents*numgrps + num_cmpds_monitored*max_cmpd_events*2 + add_elsize < eventlog_maxelements * 2
        parse_by_parts = 0;
    else    % slow parsing needed
        parse_by_parts = 1;
    end


    if run_id == 1    % only ask the question for first run.  Could flag it so this question pops up if later run (and only later run) has
                    % potential memory problem, but what if user walks away?
        if tot_elements_needed > eventlog_maxelements *2
            if parse_by_parts
                user_reply = input(['The chosen file has ' num2str(tot_events) ' reaction events.' ...
                                    '\nParsing and analysis will require memory able to handle at least ' num2str(tot_elements_needed) ' elements.' ...
                                    '\nThis exceeds the value of eventlog_maxelements*2 of ' num2str(eventlog_maxelements*2) '.' ...
                                    '\nIf you continue, out of memory errors are possible.' ...
                                    '\nThe internal files needed are also large, and will require a slower parsing method.' ...
                                    '\nContinue anyways? [y]/n: '], 's');
            else
                user_reply = input(['The chosen file has ' num2str(tot_events) ' reaction events.' ...
                                    '\nParsing and analysis will require memory able to handle at least ' num2str(tot_elements_needed) ' elements.' ...
                                    '\nThis exceeds the value of eventlog_maxelements*2 of ' num2str(eventlog_maxelements*2) '.' ...
                                    '\nIf you continue, out of memory errors are possible.  \nContinue anyways? [y]/n: '], 's');
            end
            if isempty(user_reply)
                user_reply = 'y';
            else
                user_reply = lower(user_reply);
            end
            if ~strcmp(user_reply,'y')
                clear global stopwait
                if exist('hwait')
                    delete(hwait)
                end
                out = 0;
                return
            end
        elseif parse_by_parts
            user_reply = input(['\nThe internal files needed for this parsing are large and will require a slower parsing method' ...
                                '\nin order to reduce the chance of an out-of-memory error.  If you continue to get this message with' ...
                                '\ndata sets that do not cause an out-of-memory error,' ...
                                '\nthen it might be that eventlog_maxelements is set too small in storing_and_plotting.m.' ...
                                '\nContinue anyways? [y]/n: '], 's');
            if isempty(user_reply)
                user_reply = 'y';
            else
                user_reply = lower(user_reply);
            end
            if ~strcmp(user_reply,'y')
                clear global stopwait
                if exist('hwait')
                    delete(hwait)
                end
                return
            end
        end
    end
    
    if ~parse_by_parts
        % make an oversized vector for each compound filled with zeros,
        % i.e, preallocate memory for speed
        temp_rxntimes = zeros(max_rxnum,max_rxnevents,numgrps);
        temp_data_time = zeros(num_cmpds_monitored,max_cmpd_events+1);
        temp_data_cmpds = zeros(num_cmpds_monitored,max_cmpd_events+1);
        for i = 1:num_cmpds_monitored
            n_event_cmp(i) = 1;
            temp_data_time(i,1) = time_start;
            temp_data_cmpds(i,1) = cmpd_vals(monitor_cmpds_gid(i),1);
        end
        
        for nlog = 1:num_elfiles

            ev_log_name = [evlog_part2name '_' num2str(nlog) ];
            if exist(fullfile(output_dir, ev_log_name),'file')
                event_log = read_logdata(fullfile(output_dir, ev_log_name));
            else
                error([fullfile(output_dir, ev_log_name) ' does not exist'])
            end

            num_events = size(event_log,1);


            for i = 1:num_events
                rxnnum = event_log(i,2);
                if grpnum_exists
                    grp_num = event_log(i,3);
                else
                    grp_num = 1;
                end
                

                n_event_rxn(rxnnum,grp_num) = n_event_rxn(rxnnum,grp_num) + 1;
                temp_rxntimes(rxnnum,n_event_rxn(rxnnum,grp_num),grp_num) = event_log(i,1);

                for naff = 1:num_cmpds_affected(rxnnum,grp_num)
                    this_gid = cmpds_affected(rxnnum,naff,grp_num);
                    cmpd_vals(this_gid,1) = cmpd_vals(this_gid,1) + stoich_data(rxnnum,naff,grp_num);

                    monidx = monlist_idx(rxnnum,naff,grp_num);
                    n_event_cmp(monidx) = n_event_cmp(monidx) + 1;
                    temp_data_time(monidx,n_event_cmp(monidx)) = event_log(i,1);
                    temp_data_cmpds(monidx,n_event_cmp(monidx)) = cmpd_vals(this_gid,1);
                end

            end

        end     %nlog

        for ngrp = 1:numgrps
            for nrx = 1:numrxns(ngrp)
                rxn_eventtimes{nrx,ngrp} = temp_rxntimes(nrx,1:times_rxns_executed{ngrp, run_id}(nrx),ngrp);
            end
        end

        for ncmpd = 1:num_cmpds_monitored
            data_time{ncmpd} = temp_data_time(ncmpd,1:n_event_cmp(ncmpd));
            data_cmpds{ncmpd} = temp_data_cmpds(ncmpd,1:n_event_cmp(ncmpd));
        end
   
    else    % parse by parts


        for nlog = 1:num_elfiles


            ev_log_name = [evlog_part2name '_' num2str(nlog) ];
            if exist(fullfile(output_dir, ev_log_name),'file')
                event_log = read_logdata(fullfile(output_dir, ev_log_name));
            else
                error([fullfile(output_dir, ev_log_name) ' does not exist'])
            end

            num_events = size(event_log,1);
            
            for ngrp = 1:numgrps
                for nrx = 1:numrxns(ngrp)
                    % i.e, preallocate memory for speed
                    temp_rxntimes = zeros(1,times_rxns_executed{ngrp, run_id}(nrx));
                    n_event_rxn = 0;
                    for i = 1:num_events
                        if grpnum_exists
                            grp_num = event_log(i,3);
                        else
                            grp_num = 1;
                        end
                        if grp_num == ngrp
                            rxnnum = event_log(i,2);
                            if rxnnum == nrx
                                n_event_rxn = n_event_rxn + 1;
                                temp_rxntimes(n_event_rxn) = event_log(i,1);
                            end
                        end
                    end
                    temp_rxntimes(n_event_rxn+1:end) = [];
                    rxn_eventtimes{nrx,ngrp} = [rxn_eventtimes{nrx,ngrp} temp_rxntimes];
                end
            end

            
            
            for ncmpd = 1:num_cmpds_monitored
                temp_data_time = zeros(1,each_cmpd_events(ncmpd)+1);
                temp_data_cmpds = zeros(1,each_cmpd_events(ncmpd)+1);
                this_gid = monitor_cmpds_gid(ncmpd);
                if nlog == 1
                    n_event_cmpd = 1;
                    temp_data_time(1) = time_start;
                    temp_data_cmpds(1) = cmpd_vals(this_gid,1);
                    cmpd_val = cmpd_vals(this_gid,1);
%                    tstep_iter = 1;
%                    this_montime = mon_cnums(tstep_iter).time;
                else
                    n_event_cmpd = 0;
                    cmpd_val = data_cmpds{ncmpd}(end);

                end

                
                for i = 1:num_events
                    rxnnum = event_log(i,2);
                    if grpnum_exists
                        grp_num = event_log(i,3);
                    else
                        grp_num = 1;
                    end
                    
                    
                    for naff = 1:num_cmpds_affected(rxnnum,grp_num)
                        monidx = monlist_idx(rxnnum,naff,grp_num);
                        if monidx == ncmpd
                            n_event_cmpd = n_event_cmpd + 1;
                            cmpd_val = cmpd_val + stoich_data(rxnnum,naff,grp_num);

                            temp_data_time(n_event_cmpd) = event_log(i,1);
                            temp_data_cmpds(n_event_cmpd) = cmpd_val;
                        end
                    end
                end
                temp_data_time(n_event_cmpd+1:end) = [];
                temp_data_cmpds(n_event_cmpd+1:end) = [];
                data_time{ncmpd} = [data_time{ncmpd} temp_data_time];
                data_cmpds{ncmpd} = [data_cmpds{ncmpd} temp_data_cmpds];
            end

        end     %nlog
        
    end % ~parse_by_parts

    save_mat(output_dir, output_filename, 'parsed', run_id, 0, data_time, 'data_time', data_cmpds, 'data_cmpds', rxn_eventtimes, 'rxn_eventtimes');
out = 1;
end     %i_set

clear global stopwait
if exist('hwait')
    delete(hwait)
end

else

    if exist(fullfile(output_dir, evlog_part2name),'file')
        
     for run_id = 1:num_datafiles

     evlog_name = [evlog_partname '_' num2str(run_id)];
     datafile_name =  [data_partname '_' num2str(run_id)];
    
        if exist(fullfile(output_dir, datafile_name),'file')
             [temp_data_time, temp_data_cmpds] = read_snapshot(fullfile(output_dir, datafile_name)); 
        else
             error([fullfile(output_dir, datafile_name) ' does not exist'])
        end
        
        if exist(fullfile(output_dir, evlog_name),'file')
             [temp_data_time, temp_rxntimes] = read_evlog_ag(fullfile(output_dir, evlog_name)); 
        else
             error([fullfile(output_dir, evlog_name) ' does not exist'])
        end
 
        n_event = size( temp_data_time,2);
        for ncmpd = 1:num_cmpds_monitored
            data_time{ncmpd} = temp_data_time(1,1:n_event);
            data_cmpds{ncmpd} = temp_data_cmpds(1:n_event,ncmpd);
        end
     
        ngrp = 1;
        for nrx = 1:numrxns(ngrp)
                rxn_events{nrx} = temp_rxntimes(1:n_event,nrx);
        end

        save_mat(output_dir, output_filename, 'parsed', run_id, 0, data_time, 'data_time', data_cmpds, 'data_cmpds', rxn_events,'rxn_events');
         
     end
  
    else
     error([fullfile(output_dir, evlog_part2name) ' does not exist'])   
    end
end
       