function varargout = twa_calc(varargin)

% need to implement twa.create_histplots, twa.rungrp, twa.binsize
% be sure to handle delta_t = binsize = 0
% hist_cmpds -> plotcmpds, twa.rxngrp

% used to analyze parsed eventlog files for parameter sets.


% Modified 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);

analysis1.timestart = timestart;
analysis1.timeend = timeend;


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

if ~isempty(pf)
    name_core = file_string(1:pf(1)-1);
    ssflag = 0;
    data_partname = [name_core '_parsed'];
    rundata_file = [name_core '_rundata.mat'];
elseif ~isempty(elf)
    name_core = file_string(1:elf(1)-1);
    ssflag = 0;
    data_partname = [name_core '_parsed'];
    rundata_file = [name_core '_rundata.mat'];
elseif ~isempty(df)
    name_core = file_string(1:df(1)-1);
    ssflag = 1;
    data_partname = [name_core '_data'];
    rundata_file = [name_core '_rundata.mat'];
end

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

tw = [];
cmpdnames = [];
mean_tw = [];
std_tw = [];
sens_mean = [];
sens_std = [];


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

% create list of compounds that are to be included in plots, each compound occurring only once in this list
if cmpdgrp
    cmpdlist = [];
    numfigs = length(plotcmpds);
    if isempty(plotcmpds{numfigs})
        plotcmpds(numfigs) = [];
        numfigs = numfigs-1;
    end

    % start with first figure
    cmpd_num = 1;
    cmpdlist{cmpd_num} = plotcmpds{1}(1);
    fig_idx{1} = 1;

    numcmpds = length(plotcmpds{1});
    for j = 2:numcmpds
        cmpd_num = cmpd_num + 1;
        cmpdlist{cmpd_num} = plotcmpds{1}(j);
        fig_idx{1} = [fig_idx{1} j];
    end

    % do rest of figures
    for i = 2:numfigs
        numcmpds = length(plotcmpds{i});
        for j = 1:numcmpds
            cmpdlist_idx = find(strcmp(cmpdlist,plotcmpds{i}(j)));
            if isempty(cmpdlist_idx)    % a new compound for the list
                cmpd_num = cmpd_num + 1;
                cmpdlist{cmpd_num} = plotcmpds{i}(j);
                if length(fig_idx) < i
                    fig_idx{i} = cmpd_num;
                else
                    fig_idx{i} = [fig_idx{i} cmpd_num];
                end
            else    % this compound was already added to cmpdlist
                if length(fig_idx) < i
                    fig_idx{i} = cmpdlist_idx;
                else
                    fig_idx{i} = [fig_idx{i} cmpdlist_idx];
                end
            end
        end
    end
else
    cmpdlist = plotcmpds{1};
    numfigs = length(cmpdlist);
    for i = 1:numfigs
        fig_idx{i} = i;
        plotcmpds{i} = cmpdlist(i);
    end
end

deltat = analysis1.timeend - analysis1.timestart;
if binsize == 0
    numbins = 1;
else
    numbins = ceil(deltat/binsize);
end
numcmpds_list = length(cmpdlist);

first_runid = first_psetrunid + firstrun - 1;
nruns = lastrun - firstrun + 1;

% calculate time-weighted average values for each bin of each chosen compound
analysis1.cmpd = cmpdlist;

% Run on multiple processors

i = 1;
for run_id = my_rank : comm_size : nruns-1
    myrun = first_runid + run_id;
[tw0, cmpdnames1, temp1, temp2] = twui_ave(1, myrun, data_dir, data_partname, analysis1, plot_info, ssflag, binsize);
    tw(myrun,:,:) = tw0;
    out.run(i) = myrun;
    i = i + 1;
end
    out.data = tw; 

if my_rank == master

% Collect info from the nodes

n_received = i;

        while n_received < nruns

           [message_ranks, message_tags] = MPI_Probe( '*', '*', comm);
          tw0 = [];
      % 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);
                tw(run,:,:) = tw0.data(run,:,:);
             end
                n_received = n_received + i_runs;

             end
        end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% frequency histogram part

if create_histplot
    
    tw_hist = [];
    
    for i = 1:nruns
        for k = 1:numcmpds_list
            tw_hist(i,k) = sum(tw(i,k,:));
            
            % last bin may be different size
            if binsize == 0
                adj_numbins = 1;
            else
                lastbinsize = deltat - (numbins-1)*binsize;
                frac_lastbinsize = lastbinsize / binsize;
                adj_numbins = numbins - 1 + frac_lastbinsize;
            end
            tw_hist(i,k) = tw_hist(i,k) ./ adj_numbins; % the corrected average
        end
    end
    for k = 1:numcmpds_list
        mean_tw(k) = mean(tw_hist(:,k));    % not plotted or shown, but available if somebody wants to dig it out of saved file
        std_tw(k) = std(tw_hist(:,k));
    end
end
        

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot vs. time part
if rungrp == 0 | rungrp == 1
    if nruns == 1 & numcmpds_list == 1
        bin_avg = reshape(tw, 1, numbins);
        bin_sd = zeros(1, numbins);
    elseif nruns == 1 & numcmpds_list ~= 1
        bin_avg = reshape(tw, numcmpds_list, numbins);
        bin_sd = zeros(numcmpds_list, numbins);
    else
        bin_avg = reshape(mean(tw,1),numcmpds_list,numbins);
        bin_sd = reshape(std(tw,0,1),numcmpds_list,numbins);
    end
end

if numbins == 1
    bintime(1) = analysis1.timestart + (0.5).*binsize;
else
    for i = 1:numbins-1
        bintime(i) = analysis1.timestart + (i-0.5).*binsize;
    end

    last_maxt = bintime(numbins-1) + binsize./2;
    bintime(numbins) = last_maxt + (analysis1.timeend -last_maxt)./2;
end
%bintime(numbins) = analysis1.timeend - (analysis1.timeend - bintime(numbins-1) + binsize./2)/2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


savefile = fullfile(analysis_dir, analysis_file);

clear pstruc
pstruc.plot_info = plot_info;
pstruc.run_info = run_info;
pstruc.dir_info = dir_info;
pstruc.cmpd_info = cmpd_info;
pstruc.rxn_grp = rxn_grp;
pstruc.plotcmpds = plotcmpds;
pstruc.data_partname = data_partname;
pstruc.analysis1 = analysis1;
pstruc.create_histplot = create_histplot;
pstruc.analysis_type = 'average';
pstruc.rungrp = rungrp;
pstruc.firstrun = firstrun;
pstruc.lastrun = lastrun;
pstruc.timestart = analysis1.timestart;
pstruc.timeend = analysis1.timeend;

pstruc.tw = tw;
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_histplot
    pstruc.tw_hist = tw_hist;
    pstruc.cmpdnames1 = cmpdnames1;
    pstruc.first_runid = first_runid;
    pstruc.mean_tw = mean_tw;
    pstruc.std_tw = std_tw;
    pstruc.nruns = nruns;
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');


