function varargout = bns_posthoc(varargin)

% used to analyze parsed eventlog files for parameter sets.

% If 0 arguments, the parameters directory, output_directory, and setsname are obtained from
% default_directories.m..  It will use analysis_list.m in parameters
% directory for the analysis.

% if 1 argument, the one argument must be the file name of the analysis
% file to be used.  This file must be in the parameters directory.  The
% parameters directory, output_directory, and setsname are obtained from
% default_directories.m.

% If 4 arguments, the first must be the file name of the analysis
% file to be used, the 2nd must be full path to output directory, the 3rd must be
% the setsname (currently called output_filename), and the 4th is the full path
% to the parameters directory

if nargin == 0 | nargin == 1
    default_directories
    param_dir = fullfile(models_dir,specific_model_dir,parameters_dir);
    
    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'];
    parsed_partname = [output_filename '_parsed'];
end

if nargin == 0
    analysis_file = 'analysis_list.m';
elseif nargin == 1
    analysis_file = varargin{1};
elseif nargin == 4
    analysis_file = varargin{1};
    output_dir = varargin{2};
    param_dir = varargin{4};
    data_partname = [varargin{3} '_data'];
    parsed_partname = [varargin{3} '_parsed'];
    output_filename = varargin{3};
else
    disp('Incorrect number of arguments.  Zero arguments means to use default_directories.m and analysis_list.m')
    disp('One argument is the analysis file in parameters directory.  Other information grabbed from default_directories.m')
    disp('If four arguments, they should be (1) analysis file in parameters directory, (2) full pathname of the output file, (3) output_filename, and (4) full pathname of parameters directory.')    
    return
end

[matlabversion,this_matlab_date] = version;
this_matlab_date = datenum(this_matlab_date);
matlab_700_date = datenum('July 1, 2003');
if this_matlab_date > matlab_700_date
    this_matlab_version = 7;
else
    this_matlab_version = 6;
end

curdir = pwd;
cd(param_dir);
if exist(analysis_file,'file')
    eval(analysis_file)
    cd(curdir);
else
    cd(curdir);
    error([analysis_file ' does not exist in parameters directory ' param_dir '.']);
end

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

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

if exist(fullfile(output_dir, datafile_list(1).name),'file')
    load(fullfile(output_dir, datafile_list(1).name), 'plot_info', 'run_info', 'dir_info')
else
    error([fullfile(output_dir, datafile_list(1).name) ' does not exist'])
end

pset = run_info.param_sets;

if isempty(pset)
    error('bns_posthoc only works in conjunction with parameter sets')
end

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
    

num_analyses = length(analysis);

for i = 1:num_analyses
    
    tw = [];
    cmpdnames = [];
    mean_tw = [];
    std_tw = [];
    sens_mean = [];
    sens_std = [];
    
%    first_datafile_name = datafile_list(pset(analysis(i).pset).first_runid).name;
%    if exist(fullfile(output_dir, first_datafile_name),'file')
%        load(fullfile(output_dir, first_datafile_name), 'cmpd_info', 'rxn_grp')
%    else
%        error([fullfile(output_dir, first_datafile_name) ' does not exist'])
%    end
    
    
    switch lower(analysis(i).type)
        
        case 'average'
            
            this_pset = pset(analysis(i).pset);
            num_runs = this_pset.numruns;
            first_runid = this_pset.first_runid;
            if exist(fullfile(output_dir, datafile_list(first_runid).name),'file')
                load(fullfile(output_dir, datafile_list(first_runid).name), 'cmpd_info', 'rxn_grp')
            else
                error([fullfile(output_dir, datafile_list(first_runid).name) ' does not exist'])
            end

            [tw, cmpdnames] = tw_ave(num_runs, first_runid, output_dir, parsed_partname, analysis(i), plot_info);
            for k = 1:length(cmpdnames)
                mean_tw(k) = mean(tw{k});
                std_tw(k) = std(tw{k});
            end
            
            savefile = fullfile(output_dir, [output_filename '_analysis_' num2str(i)]);
            save(savefile, 'tw', 'cmpdnames', 'mean_tw', 'std_tw', 'analysis', 'plot_info', 'run_info', 'dir_info', 'cmpd_info', 'rxn_grp');
            
            if isfield(analysis(i),'show_histograms')
                if analysis(i).show_histograms == 1;
                    for k = 1:length(cmpdnames)
                        figure
                        hist(tw{k})
                        hleg = legend(cmpdnames{k});
                        if this_matlab_version == 7;
                            set(hleg, 'Interpreter','none');
                        end
                    end
                    drawnow
                end
            end
            
            
        case 'sensitivity'
            
            for j = 1:2
                this_pset = pset(analysis(i).pset(j));
                num_runs = this_pset.numruns;
                first_runid = this_pset.first_runid;
                if exist(fullfile(output_dir, datafile_list(first_runid).name),'file')
                    load(fullfile(output_dir, datafile_list(first_runid).name), 'cmpd_info', 'rxn_grp')
                else
                    error([fullfile(output_dir, datafile_list(first_runid).name) ' does not exist'])
                end
                
                [t1, t2] = tw_ave(num_runs, first_runid, output_dir, parsed_partname, analysis(i), plot_info);
                tw{j} = t1;
                cmpdnames{j} = t2;
                for k = 1:length(cmpdnames{j})
                    mean_tw{j}(k) = mean(tw{j}{k});
                    std_tw{j}(k) = std(tw{j}{k});
                end
                
                if isfield(analysis(i),'show_histograms')
                    if analysis(i).show_histograms == 1;
                        for k = 1:length(cmpdnames{j})
                            figure
                            hist(tw{j}{k})
                            hleg = legend(cmpdnames{j}{k});
                            if this_matlab_version == 7;
                                set(hleg, 'Interpreter','none');
                            end
                        end
                        drawnow
                    end
                end
                
            end
            
            for k = 1:length(cmpdnames{1})
                num = (mean_tw{2}(k) - mean_tw{1}(k)) ./ ((mean_tw{2}(k) + mean_tw{1}(k)) ./ 2);
                den = (analysis(i).sensitive_to.val(2) - analysis(i).sensitive_to.val(1)) ./ ((analysis(i).sensitive_to.val(2) + analysis(i).sensitive_to.val(1)) ./ 2);
                sens_mean(k) = num ./ den;
                num = (std_tw{2}(k) - std_tw{1}(k)) ./ ((std_tw{2}(k) + std_tw{1}(k)) ./ 2);
                sens_std(k) = num ./ den;
            end
            
            figure
            bar(sens_mean)
            set(gca, 'xticklabel', cmpdnames{1});
            set(gca, 'FontSize', 7);
            title_line1 = ['Sensitivity of mean of each compound to ' num2str(den*100) '% variation in ' analysis(i).sensitive_to.name];
            title_line2 = ['time range = ' num2str(analysis(i).timestart) ' to ' num2str(analysis(i).timeend)];
            htitle = title({title_line1; title_line2});
            if this_matlab_version == 7;
                set(htitle, 'Interpreter','none');
            end
            
                        
            figure
            bar(sens_std)
            set(gca, 'xticklabel', cmpdnames{1});
            set(gca, 'FontSize', 7);
            title_line1 = ['Sensitivity of standard deviation of each compound to ' num2str(den*100) '% variation in ' analysis(i).sensitive_to.name];
            title_line2 = ['time range = ' num2str(analysis(i).timestart) ' to ' num2str(analysis(i).timeend)];
            htitle = title({title_line1; title_line2});
            if this_matlab_version == 7;
                set(htitle, 'Interpreter','none');
            end
            
            
            savefile = fullfile(output_dir, [output_filename '_analysis_' num2str(i)]);
            save(savefile, 'tw', 'cmpdnames', 'mean_tw', 'std_tw', 'sens_mean', 'sens_std', 'analysis', 'plot_info', 'run_info', 'dir_info', 'cmpd_info', 'rxn_grp');

                
    end
end
