%function twa_calc(twa)

% 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.


file_string = data_file;
%data_dir = twa.data_dir;
%analysis_dir = twa.analysis_dir;
%analysis_file = twa.analysis_file;
param_set_num = str2num(param_string);

%analysis1.cmpd = twa.plotcmpds;
%plotcmpds = twa.plotcmpds;
analysis1.timestart = timestart;
analysis1.timeend = timeend;

%binsize = twa.binsize;
%rungrp = twa.rungrp;
%cmpdgrp = twa.cmpdgrp;

%firstrun = twa.firstrun;
%lastrun = twa.lastrun;

%create_histplot = twa.create_histplots;

pf = findstr(file_string,'_parsed_1');
df = findstr(file_string,'_data_1');
if ~isempty(pf)
    name_core = file_string(1:pf(1)-1);
    datafile_core = [name_core '_data_*.mat'];
    ssflag = 0;
    data_partname = [name_core '_parsed'];
    rundata_file = [name_core '_rundata.mat'];
elseif ~isempty(df)
    name_core = file_string(1:df(1)-1);
    datafile_core = [name_core '_data_*.mat'];
    ssflag = 1;
    data_partname = [name_core '_data'];
    rundata_file = [name_core '_rundata.mat'];
end


datafilecore_wpath = fullfile(data_dir, datafile_core);

datafile_list = dir(datafilecore_wpath);
%num_datafiles = size(datafile_list,1);

%parsed_partname = [name_core '_parsed'];


%parsedfile_list = dir(fullfile(output_dir, [parsed_partname '*.mat']));
%num_parsedfiles = size(parsedfile_list,1);

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


adatafilename_thispset = [name_core '_data_' num2str(first_psetrunid) '.mat'];

% 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;
[tw, cmpdnames1, temp1, temp2] = twui_ave(nruns, first_runid, data_dir, data_partname, analysis1, plot_info, ssflag, binsize);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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')

if use_batch == 0
    plot_twa(pstruc)
end



