function varargout = reaction_freq_analysis(varargin)

% used to analyze parsed eventlog files for parameter sets.
% Created to work in MPI mode v3.1
%                            Slava 09/07   

tic
% Initialize MPI
MPI_Init;

global MPI_COMM_WORLD

% Create communicator
comm = MPI_COMM_WORLD;

% Get size and rank
comm_size = MPI_Comm_size(comm);
my_rank = MPI_Comm_rank(comm);

%Set the master node
master = 0;
% Create base message tegs
 coefs_tag = 10000;
 input_tag = 20000;
 output_tag = 30000;

analysis_constants


file_string = data_file;
param_set_num = str2num(param_string);
plotrxns = plotcmpds;

pf = findstr(file_string,'_parsed_1');
df = findstr(file_string,'_data_1');
elf = findstr(file_string,'_eventlog_1');

if ~isempty(pf)
    name_core = file_string(1:pf(1)-1);
elseif ~isempty(elf)
    name_core = file_string(1:elf(1)-1);
elseif ~isempty(df)
    name_core = file_string(1:df(1)-1);
end

    rundata_file =  [name_core '_rundata.mat'];
    parsed_partname = [name_core '_parsed'];
    data_partname = [name_core '_data'];

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


pset = run_info.param_sets;

if ~isempty(pset)
    num_psets = length(pset);
    cur_runid = 1;
    for i = 1:num_psets
        pset(i).first_runid = cur_runid;
        cur_runid = cur_runid + pset(i).numruns;
    end
    this_pset = pset(param_set_num);
    num_runs = this_pset.numruns;
    first_psetrunid = this_pset.first_runid;
else
    num_runs = run_info.num_runs;
    first_psetrunid = 1;
end


num_rxns = [];
first_runid = first_psetrunid + firstrun - 1;
nruns = lastrun - firstrun + 1;
cur_runid = first_runid;
num_grps = length(rxn_grp);
for i = 1:num_grps
    num_rxns(i) = length(rxn_grp{i}.rxn);
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if create_barchart
    
    for i = 1:num_grps
        for j = 1:num_rxns(i)
            eventnum_matrix(i).rxn(j).num = [];
        end
    end
    
    for k = 1:nruns
% Modified to work with snapshot data
%                Slava  08/07
       if ~isempty(pf)
        this_filename = [parsed_partname '_' num2str(cur_runid) '.mat'];
        this_file_wpath = fullfile(data_dir, this_filename);
        if exist(this_file_wpath,'file')
            load(this_file_wpath)
        else
            error([this_file_wpath ' does not exist'])
        end
        agregated = 0; 

            for i = 1:num_grps
                for j = 1:num_rxns(i)
                    tidx = find(rxn_eventtimes{j,i} >= timestart & rxn_eventtimes{j,i} <= timeend);
                    numevents = length(tidx);
                    eventnum_matrix(i).rxn(j).num = [eventnum_matrix(i).rxn(j).num numevents];
                end
            end
        else
% Work with agregated number of reaction
%                         Slava   08/07
        this_filename = [data_partname '_' num2str(cur_runid)];
        this_file_wpath = fullfile(data_dir, this_filename);
        if exist(this_file_wpath,'file')
            [data_time, data_cmpds, rxn_events] = read_snapshot(this_file_wpath);  
        else
            error([this_file_wpath ' does not exist'])
        end
        agregated = 1; 
            for i = 1:num_grps
                for j = 1:num_rxns(i)
                    tidx = rxn_events(:,j);
                    numevents = sum(tidx);
                    eventnum_matrix(i).rxn(j).num = [eventnum_matrix(i).rxn(j).num numevents];
                end
            end
        end
        cur_runid = cur_runid + 1;        
    end
    
    cur_rxn_num = 0;
    for i = 1:num_grps
        for j = 1:num_rxns(i)
            cur_rxn_num = cur_rxn_num + 1;
            avg_rnum(cur_rxn_num) = mean(eventnum_matrix(i).rxn(j).num);
            std_rnum(cur_rxn_num) = std(eventnum_matrix(i).rxn(j).num);
            rxn_names{cur_rxn_num} = rxn_grp{i}.rxn{j}.name;
        end
    end
    
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% reaction counts vs. time
if length(plotrxns) ~= 0 && ~isempty(plotrxns{1})
    % analysis section

    % create list of reactions that are to be included in plots, each reaction occurring only once in this rxnlist
    if rxngrp
        rxnlist = [];
        numfigs = length(plotrxns);
        if isempty(plotrxns{numfigs})
            numfigs = numfigs-1;
        end

        % start with first figure
        rxn_num = 1;
        rxnlist{rxn_num} = plotrxns{1}(1);
        fig_idx{1} = 1;

        numrxns_thisfig = length(plotrxns{1});
        for j = 2:numrxns_thisfig
            rxn_num = rxn_num + 1;
            rxnlist{rxn_num} = plotrxns{1}(j);
            fig_idx{1} = [fig_idx{1} j];
        end

        % do rest of figures
        for i = 2:numfigs
            numrxns_thisfig = length(plotrxns{i});
            for j = 1:numrxns_thisfig
                rxnlist_idx = find(strcmp(rxnlist,plotrxns{i}(j)));
                if isempty(rxnlist_idx)    % a new compound for the list
                    rxn_num = rxn_num + 1;
                    rxnlist{rxn_num} = plotrxns{i}(j);
                    if length(fig_idx) < i
                        fig_idx{i} = rxn_num;
                    else
                        fig_idx{i} = [fig_idx{i} rxn_num];
                    end
                else    % this compound was already added to cmpdlist
                    if length(fig_idx) < i
                        fig_idx{i} = rxnlist_idx;
                    else
                        fig_idx{i} = [fig_idx{i} rxnlist_idx];
                    end
                end
            end
        end
    else
        rxnlist = plotrxns{1};
        numfigs = length(rxnlist);
        for i = 1:numfigs
            fig_idx{i} = i;
            plotrxns{i} = rxnlist(i);
        end
    end
        
    
    deltat = timeend - timestart;
    numbins = ceil(deltat/binsize);
    numrxns_list = length(rxnlist);
    
    % find sim reaction numbers and rxn_grp for the reactions chosen in gui
    r_index = [];
    for m = 1:numrxns_list
        idx_found = 0;
        for i = 1:num_grps
            for j = 1:num_rxns(i)
                if strcmp(rxnlist{m}, rxn_grp{i}.rxn{j}.name) == 1
                   r_index(m,:) = [i j];
                   idx_found = 1;
                   break
                end
            end
            if idx_found
                break
            end
        end
    end
    
%    clear evmat
    binnum = zeros(nruns,numrxns_list,numbins);
    bintime = zeros(numbins,1);
    clear bintimes

    % Run on multiple processors

k = 1;
for run_id = my_rank : comm_size : nruns-1
     myrun = first_runid + run_id;
    
      if ~isempty(pf)
        this_filename = [parsed_partname '_' num2str(myrun) '.mat'];
        this_file_wpath = fullfile(data_dir, this_filename);
        if exist(this_file_wpath,'file')
            load(this_file_wpath)
        else
            error([this_file_wpath ' does not exist'])
        end

        agregated = 0; 
          for m = 1:numrxns_list
            rgrp = r_index(m,1);
            rnum = r_index(m,2);
            
            tidx = find(rxn_eventtimes{rnum,rgrp} >= timestart & rxn_eventtimes{rnum,rgrp} <= timeend);
            numevents = length(tidx);
            this_evtimes = rxn_eventtimes{rnum,rgrp}(tidx);

            if my_rank == master
                if create_eventplots & k == 1
                    bintimes(m).times = this_evtimes;
                end
            end
            
            cur_maxt = timestart + binsize;
            if cur_maxt > timeend
                cur_maxt = timeend
            end
            prev_maxt = timestart;
            curbin = 1;
            bincount = 0;
            n = 1;
            while n <= numevents
                curtime = this_evtimes(n);
                if curtime <= cur_maxt
                    bincount = bincount + 1;
                    n = n + 1; 
                end
                if curtime > cur_maxt | n > numevents
                    if cur_maxt == timeend
                        bin_deltat = timeend - prev_maxt;
                    else
                        bin_deltat = binsize;
                    end
                    bintime(curbin) = cur_maxt - bin_deltat/2;
                    binnum(run_id,m,curbin) = bincount./bin_deltat;
                    bincount = 0;
                    if curtime <= timeend
                        curbin = curbin + 1;
                        prev_maxt = cur_maxt;
                        cur_maxt = cur_maxt + binsize;
                        if cur_maxt > timeend
                            cur_maxt = timeend;
                        end
                    end
                end % curtime
            end  % while 
          end %rxnlist
    
        else

 % Work with agregated number of reaction
%                         Slava   08/07
        this_filename = [data_partname '_' num2str(myrun)];
        this_file_wpath = fullfile(data_dir, this_filename);
        if exist(this_file_wpath,'file')
            [data_time, data_cmpds, rxn_events] = read_snapshot(this_file_wpath);  
        else
            error([this_file_wpath ' does not exist'])
        end
        
        agregated = 1; 
        
          for m = 1:numrxns_list
            rgrp = r_index(m,1);
            rnum = r_index(m,2);
            
            tidx = rxn_events(:,rnum);
            numevents = length(tidx);
            this_evtimes = data_time(1,1);
            
            if my_rank == master
                if create_eventplots & k == 1
                    bintimes(m).times = data_time;
                    bintimes(m).agnum = tidx;
                end
            end

            snapshot_interval =  data_time(1,2) - data_time(1,1);
            if (binsize < snapshot_interval)
                error( ' Time binsize is smaller then snapshot interval!!! ')
            end
        
            cur_maxt = timestart + binsize;
            if cur_maxt > timeend
                cur_maxt = timeend;
            end
            prev_maxt = timestart;
            curbin = 1;
            bincount = 0;
            n = 1;
            while n <= numevents
                curtime = data_time(1,n);
                if curtime <= cur_maxt
                    bincount = bincount + tidx(n);
                    n = n + 1; 
                end
                if curtime > cur_maxt | n > numevents
                    if cur_maxt == timeend
                        bin_deltat = timeend - prev_maxt;
                    else
                        bin_deltat = binsize;
                    end
                    bintime(curbin) = cur_maxt - bin_deltat/2;
                    binnum(run_id,m,curbin) = bincount./bin_deltat;
                    bincount = 0;
                    if curtime <= timeend
                        curbin = curbin + 1;
                        prev_maxt = cur_maxt;
                        cur_maxt = cur_maxt + binsize;
                        if cur_maxt > timeend
                            cur_maxt = timeend;
                        end
                    end
                end % curtime
            end  % while          
          end %rxnlist
        end  % if/else 
     out.run(k) = myrun;
     k = k+1;
    end % run_id = 1:nruns

end
    out.data = binnum;  
    

if my_rank == master

% Collect info from the nodes

n_received = k;

        while n_received < nruns

           [message_ranks, message_tags] = MPI_Probe( '*', '*', comm);

      % If message_ranks is not empty - receive the output
             if ~isempty(message_ranks)

      % Sort the received messages
                [m_tags_sorted, m_tags_sorted_idx] = sort(message_tags);

                dest = message_ranks(m_tags_sorted_idx(1));
                master_tag = message_tags(m_tags_sorted_idx(1));
                mesg_num = master_tag - output_tag;

                tw0 = MPI_Recv( dest, master_tag, comm);

                i_runs = size(tw0.run,2)

             for i=1:i_runs
                run = tw0.run(i)
                binnum(run,:,:) = tw0.data(run,:,:);
             end

                n_received = n_received + i_runs;

             end
        end
    
    if rungrp == 0 | rungrp == 1
        if nruns == 1 & numrxns_list == 1
            bin_avg = reshape(binnum, 1, numbins);
            bin_sd = zeros(1, numbins);
        elseif nruns == 1 & numrxns_list ~= 1
            bin_avg = reshape(binnum, numrxns_list, numbins);
            bin_sd = zeros(numrxns_list, numbins);
        else
            bin_avg = reshape(mean(binnum),numrxns_list,numbins);
            bin_sd = reshape(std(binnum),numrxns_list,numbins);
        end
    end

    
end

savefile = fullfile(analysis_dir, analysis_file);

pstruc.plot_info = plot_info;
pstruc.run_info = run_info;
pstruc.dir_info = dir_info;
pstruc.rxn_grp = rxn_grp;
pstruc.create_barchart = create_barchart;
pstruc.plotrxns = plotrxns;
pstruc.rungrp = rungrp;
pstruc.firstrun = firstrun;
pstruc.lastrun = lastrun;
pstruc.timestart = timestart;
pstruc.timeend = timeend;
pstruc.parsed_partname = parsed_partname;
pstruc.data_partname = data_partname;
pstruc.create_eventplots = create_eventplots;
pstruc.analysis_type = 'reaction frequency';
pstruc.agregated = agregated;

if length(plotrxns) ~= 0 && ~isempty(plotrxns{1})
    
    pstruc.binnum = binnum;
    pstruc.bintime = bintime;
    pstruc.numfigs = numfigs;
    pstruc.fig_idx = fig_idx;
    pstruc.binsize = binsize;
    if rungrp == 0 | rungrp == 1
        pstruc.bin_avg = bin_avg;
        pstruc.bin_sd = bin_sd;
    end

    if create_eventplots
        pstruc.numrxns_list = numrxns_list;
        pstruc.bintimes = bintimes;
    end

end
if create_barchart
    pstruc.avg_rnum = avg_rnum;
    pstruc.std_rnum = std_rnum;
    pstruc.rxn_names = rxn_names;
    pstruc.eventnum_matrix = eventnum_matrix;
end

save(savefile,'pstruc');

else
%  Send tw-array to the master node

      master_tag = output_tag + my_rank;
      MPI_Send(master, master_tag, comm, out);

end

MPI_Finalize;
disp(' SUCCESS');




