function batch_analyse_stngpe(e_fname,flags_fname,pathroot,an_fname)

% BATCH_ANALYSE_STNGPE to analyse firing statistics of sampled STN and GPe neurons
%   BATCH_ANALYSE_STNGPE(E,F,A) analyses the extracted data in file E (a
%   string) using the analysis methods specifed in flags file F (a string).
%   Saves the results to the file named A.
%
% Computes in order for STN and GPe:
% (0) Plots raster and single-neuron trace, if requested
% (1) Mean, SD, SE, and CV of firing rates over whole simulation (includes SNr)
% (2) ISIs, ISI histograms (including input)
% (3) Burst firing analysis (KBSTA method)
% (4) Power-spectra (multi-taper method - chronux toolbox)
%     (a) finds most significant peak
%     (b) finds any significant peak in LFO (<1.5Hz) range (Magill et al. 2001)
% (5) Intra- and inter-nucleus cross-correlation, using peak/trough detection method set at start
% (6) Auto-correlograms
% (7) Coherency (chronux toolbox)
% (8) Spike-triggered averaging (including computation of pseudo-EEG on which it's based)
%
% Mark Humphries & Kevin Gurney 6/1/2006. 

load([pathroot e_fname])       % load extracted results

eval(flags_fname)   % evaluate analysis flag file

 [a s] =  system('hostname');
 host = findstr(s, 'iceberg');
 if ~isempty(host) % on iceberg
     ice_path1 = genpath('/home1/pc/pc1mdh/BG spiking model');
     ice_path2 = genpath('/home1/pc/pc1mdh/Matlab Tools');
     path(path, ice_path1);
     path(path, ice_path2);
 end

cl_fig

do_display = 0; % set for showing graphs interactively

% t_out_t = double(out_t) .* dt; % timestamps for network events
t_in_t = double(in_t) .* dt;    % timestamps for input streams
time_steps = time_seconds / dt;

clear link_n link_wAMPA link_wNMDA link_wGABAa

%%%%%%%%%%%%%%%% SET ANALYSIS PARAMETERS HERE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% parameters for ISI data
ISIbins = 100; %No. of bins in ISI hist.
limits = [0.005 0.6];  % ISI hist limits in seconds (outside these bounds get put into end bins)

% parameters for burst analysis
eta_on = 0.4;
eta_off = 0.4;

% parameters for periodogram
freq_range = [0.1 100];    % min and max freqs in periodogram
tapers = [3 5];                         % recommended in Chronuz; 5 tapers with bandwidth*time = 3
pad = 2;                                % Padding the frequency sampling by another 2 powers of 2;
Fs = 1 ./ dt;                           % best to use the underlying sampling;
err_bars = [1 0.01];                    % method 1 (theoretical) and p_sig = 0.01
fscorr = 0;                             %don't use finite size corrections
Min_spikes = 12;                        % have to have at least this number of spikes to do spectrum

LFO_freq = 1.5; % (For Magill's analysis)

% parameters for auto- and cross-correlogram
binsize = 0.001;   % in seconds
maxlag = 0.5;       % in seconds
acf_maxlag = 2;     % window size for auto-correlograms

p = 0.01;           % 99% confidence interval
num_bins = 2*maxlag/binsize + 1;
method1 = 0; % method 1 (proportion of bins outside confidence limits) or method 2 (Raz et al 2000)
if method1
    num_sig_bins = round(num_bins * p); % number of bins over confidence interval required - 
    % use for testing any bins
else
    num_sig_bins = 3;           % when using consecutive bins, more than 3 is significant?                                    
end

% parameters for coherence
freq_rangeC = [0.1 100];    % min and max freqs in periodogram
tapersC = [3 5];                         % recommended in Chronuz; 5 tapers with bandwidth*time = 3
padC = 2;                                % Padding the frequency sampling by another 2 powers of 2;
FsC = 1 ./ dt;                           % best to use the underlying sampling;
err_barsC = [2 0.01];                    % method (jacknife) and p_sig = 0.01
fscorrC = 0;                             %don't use finite size corrections
Min_spikesC = 12;                        % have to have at least this number of spikes to do spectrum


%%% parameters for spike-triggered averaging
    
% a. parameters for smoothed input firing rate estimation
    smooth_win_size = 0.01;     % width in seconds
    step_size = 10;             % window step in time-steps
    window = 'alpha';           % alt: 'gaussian'
% b. parameters for spike-triggered average
    trig_win_half = 1;          % half-size of averaging window in seconds

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

%%%%%%% some basic pre-processing %%%%%%
input_times = cell(neurons_per_nucleus,1);

% create time-stamp arrays for each input
for loop=1:neurons_per_nucleus    
    input_times{loop} = dt .* double(in_t(in_n == EXT(loop)));
end

%%%% create data structures for periodogram results - remove this later if
%%%% we tidy up data storage!

% number of neurons from each structure
n_STN = length(STN_times);
n_GPe = length(GPe_times);
n_GPi = nan;

for loop = 1:n_STN
    STN_data(loop) = struct('times', STN_times{loop}');     % for multi_taper - if multichannel use is made of the
                                                            % chronux call, then need this format for data
end

for loop = 1:n_GPe
    GPe_data(loop) = struct('times', GPe_times{loop}');   % for multi_taper
end

%%%%%%%%%% RASTER PLOT %%%%%%%%%%%%%%%%%%
if do_display & do_raster
    structured_raster(in_n, in_t, out_n, out_t, dt, neurons_per_nucleus, SD1, SD2, STN, GPe, GPi, time_seconds);
    
    display('press any key to continue');
    pause;
end

%%%%%%%%%%%% trace %%%%%%%%%%%%%%%%%%%
% show trace if appropriate
if any([trace_n == STN trace_n == GPe]) & do_trace & do_display
	trace_times = 1:1:time_steps;
	figure
	plot(trace_times,trace_vals(:,1))
	figure
	plot(trace_times,trace_vals(:,2))
	figure
	plot(trace_times,trace_vals(:,3))
    tile
    drawnow
end

%%%%%%%% basic mean firing rate stats
if do_means
    
    % determine mean outputs (firing rates) for each neuron
    mean_STN = cellfun('length',STN_times) ./ time_seconds;     % length of each array is number of spikes! 
    mean_GPe = cellfun('length',GPe_times) ./ time_seconds;
   
    
    % find mean across neurons in each nucleus
    STN_Hz = mean(mean_STN);
    GPe_Hz = mean(mean_GPe);
    
    % find std errors of firing rates in each nucleus
    std_STN = std(mean_STN);
    std_GPe = std(mean_GPe);
    
    
    sem_STN = std_STN / sqrt(neurons_per_nucleus);
    sem_GPe = std_GPe / sqrt(neurons_per_nucleus);
    
    % find CV of firing rates for each nucleus (*not* neuron)
    % for analysis of single neuron trains, find CV of ISIs...
    CV_STN = std_STN ./ STN_Hz;
    CV_GPe = std_GPe ./ GPe_Hz;
    
    if exist('GPi_times')
        n_GPi = length(GPi_times);
        mean_GPi = cellfun('length',GPi_times) ./ time_seconds;
        GPi_Hz = mean(mean_GPi);         % included for tonic firing rate setting
        std_GPi = std(mean_GPi);
        sem_GPi = std_GPi / sqrt(neurons_per_nucleus) ;       % computes the standard error of the mean (SEM) as quoted in e.g Urbain et al 2000
        CV_GPi = std_GPi ./ GPi_Hz;
    else
        n_GPi = nan;
        mean_GPi = nan;
        GPi_Hz = nan;
        std_GPi = nan;
        sem_GPi = nan;
        CV_GPi = nan;
    end 
    
%     fprintf(1, 'mean STN rate %.2f:  stderr of mean STN rate %.2f: std dev of STN rate %.2f\n', STN_Hz, sem_STN, std_STN);
%     fprintf(1, 'mean GPe rate %.2f:  stderr of mean GPe rate %.2f: std dev of STN rate %.2f\n', GPe_Hz, sem_GPe, std_GPe);
%     fprintf(1, 'mean GPi rate %.2f:  stderr of mean GPi rate %.2f: std dev of STN rate %.2f\n', GPi_Hz, sem_GPi, std_GPi);

end
if do_display
    display('press any key to continue 2');
    pause;
end
 


%%%%%%%%%%%%%%%% isi %%%%%%%%%%%%%%%%%%%%%
STN_rates = cell(n_STN,1);            % instantaneous rates for *all* events (not just bursts)
STN_x_isi = cell(n_STN,1);            % time stamps for rates above (almost just timestamps )
STN_ISIhist = cell(n_STN,1);          % ISI histograms

GPe_rates = cell(n_GPe,1);
GPe_x_isi = cell(n_GPe,1);
GPe_ISIhist = cell(n_GPe,1);

input_rates = cell(length(EXT),1);
input_x_isi = cell(length(EXT),1);
input_ISIhist = cell(length(EXT),1);

cum_STN_spikes =  sparse(1,time_steps);
cum_GPe_spikes =  sparse(1,time_steps);
cum_input_spikes =  sparse(1,time_steps);
if do_isi
    for loop = 1:n_STN
        % generate ISI values
        [STN_rates{loop},STN_ISIhist{loop},STN_x_isi{loop},STN_x_hist] = LIF_ISI_analysis(STN_times{loop},ISIbins,limits);
        STN_discrete_times = sparse(1,time_steps);
        STN_discrete_times(round(STN_times{loop} ./ dt)) = 1;
        cum_STN_spikes = cum_STN_spikes + STN_discrete_times;
    end

    for loop = 1:n_GPe
        [GPe_rates{loop},GPe_ISIhist{loop},GPe_x_isi{loop},GPe_x_hist] = LIF_ISI_analysis(GPe_times{loop},ISIbins,limits);
        GPe_discrete_times = sparse(1,time_steps);
        GPe_discrete_times(round(GPe_times{loop} ./ dt)) = 1;
        cum_GPe_spikes = cum_GPe_spikes + GPe_discrete_times;
    end
    
    for loop = 1:neurons_per_nucleus
        [input_rates{loop},input_ISIhist{loop},input_x_isi{loop},input_x_hist] = LIF_ISI_analysis(input_times{loop},ISIbins,limits);
        input_discrete_times = sparse(1,time_steps);
        input_discrete_times(round(input_times{loop} ./ dt) + 1) = 1; % input time can be zero causes a problem
        input_discrete_times = input_discrete_times(1:time_steps);
        cum_input_spikes = cum_input_spikes + input_discrete_times;
    end
    rs_STN = find_mean_isis(cum_STN_spikes);
    rs_GPe = find_mean_isis(cum_GPe_spikes);
    rs_in = find_mean_isis(cum_input_spikes);
    
    
    %%%% show individual isi pots for STN %%%%%
    No_subplots_isi = 4;
    if do_display
        done = 0;
        loop = 0;
        while loop < n_STN
            plot_number = mod(loop, No_subplots_isi) + 1;
            if  plot_number == 1
                figure
            end
            subplot(No_subplots_isi, 1, plot_number);
            plot(STN_x_isi{loop+1}, STN_rates{loop+1});
            drawnow;
            loop = loop + 1;
        end
    end
    
    if do_display
       display('press any key to continue');
       pause
    end
    %%%%%%%%%%%%%%%%%%%%
    
    
    
    if do_display
        figure
        plot(rs_STN);
        title('mean isis for STN')
        figure
        plot(rs_in);
        title('mean isis for input')
    end
end



%%%%%%%%%%%%%% bursts %%%%%%%%%%%%%%%
if do_bursts
    % determine burst status
    STN_bursts = cell(n_STN,1);           % burst start and end times
    STN_burst_isis = cell(n_STN,1);       % within burst ISIs (entire sequence for each burst)
    
    
    % repeat for GPe....
    GPe_bursts = cell(n_GPe,1);           % burst start and end times
    GPe_burst_isis = cell(n_GPe,1);       % within burst ISIs
    
    
    
    for loop = 1:n_STN
        % burst-detection
        fprintf(1, 'STN burst loop count %d\n', loop);
        
        [STN_bursts{loop},STN_burst_isis{loop}] = kbsta(STN_times{loop},[],eta_on,eta_off);
        [N_STN_bursts(loop) c] = size(STN_bursts{loop});
    end
    
    for loop = 1:n_GPe
        % burst-detection
        fprintf(1, 'GPe burst loop count %d\n', loop);
       
        [GPe_bursts{loop},GPe_burst_isis{loop}] = kbsta(GPe_times{loop},[],eta_on,eta_off);
        [N_GPe_bursts(loop) c] = size(GPe_bursts{loop});        
    end
    
    if do_display
        figure
        subplot(2,1,1)
        hist(N_STN_bursts, [0:15]);
        title('histogram of STN bursts');
        ylabel('No of cells');
        xlabel('No bursts');
        
        subplot(2,1,2)
        hist(N_GPe_bursts, [0:15]);
        title('histogram of GPe bursts');
        ylabel('No of cells');
        xlabel('No bursts');
        
        drawnow
    end

end

%%%%%%%%%%% spectra %%%%%%%%%%%%%%%%%%

% periodogram analysis
% [STN_results,STN_store_power,STN_fs_array] = scargle_analysis(STN_rates,STN_x_isi,dt, time_seconds, [], freqs);
% STNnans = isnan(STN_results(:,1)) | STN_results(:,1)==0;
% STN_sum_results = STN_results;
% STN_sum_results(STNnans,:) = [];       % results matrix without empty records
% 
% [GPe_results,GPe_store_power,GPe_fs_array] = scargle_analysis(GPe_rates,GPe_x_isi,dt, time_seconds, [],freqs);
% GPenans = isnan(GPe_results(:,1)) | GPe_results(:,1)==0;
% GPe_sum_results = GPe_results;
% GPe_sum_results(GPenans,:) = [];       % results matrix without empty records

if do_display
    display('press any key to continue 3');
    pause;
end
if do_spectra
    % multi-taper analysis
    
    STN_spect_res = struct('freqs', [], 'powers', [], 'errs', [], 'peak_freq', [], 'peak_power', [], 'LFO', []);
    GPe_spect_res = struct('freqs', [], 'powers', [], 'errs', [], 'peak_freq', [], 'peak_power', [], 'LFO', []);
    
    disp('doing multi-taper spectral analysis');
    
    for loop = 1:n_STN
        fprintf(1, 'STN spectral loop count %d\n', loop);
        
        %% STN%%% 
        data = STN_data(loop).times;
        [r c] = size(data);
        if r >= Min_spikes
            [STN_MT_power, STN_MT_fs_array, R, STN_MT_errs] = mtspectrumpt(data, tapers, pad, Fs, freq_range, err_bars, 0, fscorr); % 0 is no trial average
            STN_err_lo = STN_MT_errs(1,:);
            STN_spect_res(loop).freqs = STN_MT_fs_array;
            STN_spect_res(loop).powers = STN_MT_power;
            STN_spect_res(loop).errs = STN_MT_errs;
            
            sig_test_STN_power = STN_err_lo - mean(STN_MT_power);
            [STN_peak STN_peak_index] = max(sig_test_STN_power);
            
            if STN_peak >0
                STN_spect_res(loop).peak_freq = STN_MT_fs_array(STN_peak_index);
                STN_spect_res(loop).peak_power = STN_peak;
            else
                STN_spect_res(loop).peak_freq = NaN;
                STN_spect_res(loop).peak_power = NaN;
            end

            %%%% extra for LFOs in magill work %%%%%%%%%

            LFO_freq_indices = find(STN_MT_fs_array < LFO_freq);
            STN_LFO_sig_test = sig_test_STN_power(LFO_freq_indices);
            [STN_peak_LFO STN_peak_index_LFO] = max(STN_LFO_sig_test);
            if STN_peak_LFO > 0
                STN_spect_res(loop).LFO = 1;
            else
                STN_spect_res(loop).LFO = 0;
            end

	    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        else
            STN_spect_res(loop).freqs = [];
            STN_spect_res(loop).powers = [];
            STN_spect_res(loop).errs = [];
            STN_spect_res(loop).peak_freq = NaN;
            STN_spect_res(loop).peak_power = NaN;
            STN_spect_res(loop).LFO = NaN;
        end
    end
    
    for loop = 1:n_GPe
        fprintf(1, 'GPe spectral loop count %d\n', loop);
        %% GPe %%% 
        data = GPe_data(loop).times;
        [r c] = size(data);
        if r >= Min_spikes
            [GPe_MT_power, GPe_MT_fs_array, R, GPe_MT_errs] = mtspectrumpt(data, tapers, pad, Fs, freq_range, err_bars, 0, fscorr); % 0 is no trial average
            GPe_err_lo = GPe_MT_errs(1,:);
            GPe_spect_res(loop).freqs = GPe_MT_fs_array;
            GPe_spect_res(loop).powers = GPe_MT_power;
            GPe_spect_res(loop).errs = GPe_MT_errs;
            
            sig_test_GPe_power = GPe_err_lo - mean(GPe_MT_power);
            [GPe_peak GPe_peak_index] = max(sig_test_GPe_power);
            
            if GPe_peak >0
                GPe_spect_res(loop).peak_freq = GPe_MT_fs_array(GPe_peak_index);
                GPe_spect_res(loop).peak_power = GPe_peak;
            else
                GPe_spect_res(loop).peak_freq = NaN;
                GPe_spect_res(loop).peak_power = NaN;
            end

            %%%% extra for LFOs in magill work %%%%%%%%%

            LFO_freq_indices = find(GPe_MT_fs_array < LFO_freq);
            GPe_LFO_sig_test = sig_test_GPe_power(LFO_freq_indices);
            [GPe_peak_LFO Gpe_peak_index_LFO] = max(GPe_LFO_sig_test);
            if GPe_peak_LFO > 0
                GPe_spect_res(loop).LFO = 1;
            else
                GPe_spect_res(loop).LFO = 0;
            end

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

        else
            GPe_spect_res(loop).freqs = [];
            GPe_spect_res(loop).powers = [];
            GPe_spect_res(loop).errs = [];
            GPe_spect_res(loop).peak_freq = NaN;
            GPe_spect_res(loop).peak_power = NaN;
            GPe_spect_res(loop).LFO = NaN;
        end
    end
    
    sig_STN_peak_freqs = [STN_spect_res.peak_freq]; 
    STN_sig_fs = sig_STN_peak_freqs(find(sig_STN_peak_freqs > 0));
    lows = STN_sig_fs(STN_sig_fs < 5);
    sig_GPe_peak_freqs = [GPe_spect_res.peak_freq];
    GPe_sig_fs = sig_GPe_peak_freqs(find(sig_GPe_peak_freqs > 0));
    lowg = GPe_sig_fs(GPe_sig_fs < 5);

    if do_display
        figure
       
        subplot(4,1,1)
        hist(STN_sig_fs, 20);
        title('significant maximum  frequencies - STN');
        xlabel('peak freq');
        subplot(4,1,2)
        hist(lows, 15);
        title('significant maximum  frequencies in range 0-5Hz - STN');
        xlabel('peak freq');
        
        subplot(4,1,3)
        hist(GPe_sig_fs, 20);
        title('significant maximum  frequencies - Gpe');
        xlabel('peak freq');
        subplot(4,1,4)
        hist(lowg, 15);
        title('significant maximum  frequencies in range 0-5Hz - GPe');
        xlabel('peak freq');
        
        drawnow
    end    
end % do_spectra

%%%%%%%%%%%%%% auto- and cross-correlograms  %%%%%%%%%%
n_STN_pairs = n_STN ^ 2;
n_GPe_pairs = n_GPe ^ 2;
n_inter_pairs = n_STN * n_GPe;

STN_pairs = list_of_pairs(1:n_STN,1:n_STN,'su'); 
GPe_pairs = list_of_pairs(1:n_GPe,1:n_GPe,'su');
STNGPe_pairs = list_of_pairs(1:n_STN,1:n_GPe,'su');

STN_ccf = cell(n_STN_pairs,1);
GPe_ccf = cell(n_GPe_pairs,1);
STNGPe_ccf = cell(n_inter_pairs,1);
total_STN_ccf = zeros(1,num_bins);
total_GPe_ccf = zeros(1,num_bins);
total_STNGPe_ccf = zeros(1,num_bins);
total_STN_pairs = 0;
total_GPe_pairs = 0;
total_STNGPe_pairs = 0;

pos_sigSTNpairs = zeros(n_STN_pairs,1);
neg_sigSTNpairs = zeros(n_STN_pairs,1);
pos_sigGPepairs = zeros(n_GPe_pairs,1);
neg_sigGPepairs = zeros(n_GPe_pairs,1);


if do_xcorr_intra
disp('doing xcorr - intra nucleus');
    % OPTION 2: calculate cross-correlation histograms direct from spike trains

    GPe_peak_and_trough = 0; STN_peak_and_trough = 0; 
    for i = 1:n_STN_pairs
        fprintf(1, 'STN xcorr loop %d\n', i);

        [STN_ccf{i},xSTN,f1,f2,n_pairs] = LIF_xcorr(STN_times{STN_pairs(i,1)},STN_times{STN_pairs(i,2)},binsize,[0 time_seconds],maxlag); 
        total_STN_pairs = total_STN_pairs + n_pairs;
        [mSTN,sSTN,spSTN,stSTN] = xcorr_stats(STN_ccf{i},n_pairs,p);
 
        % test for significance in cross-correlograms
        S_pos_sig_bins = find(STN_ccf{i} >= spSTN);
        S_neg_sig_bins = find(STN_ccf{i} < stSTN);
        
        if method1
            % METHOD 1: number of bins greater than confidence interval exceeds the expected amount
            if length(S_pos_sig_bins) > num_sig_bins     % for just comparing all bins over correlogram
                pos_sigSTNpairs(i) = 1;    
            end
            if length(S_neg_sig_bins) > num_sig_bins 
                neg_sigSTNpairs(i) = 1;   
            end            
        else
            % METHOD 2: number of consecutive significant bins is greater than fixed threshold (3 - Raz et al 2000)
            for loop2 = 1:num_sig_bins
                S_pos_sig_bins = find(diff(S_pos_sig_bins) == 1); % indices of adjacent bins - reduce                
                S_neg_sig_bins = find(diff(S_neg_sig_bins) == 1); % indices of adjacent bins - reduce      
            end
            
            if ~isempty(S_pos_sig_bins)  pos_sigSTNpairs(i) = 1; end   % then there is at least one occurence of a num_sig_bins number of consecutive significant bins
            if ~isempty(S_neg_sig_bins)  neg_sigSTNpairs(i) = 1; end
        end
       
        % omit correlograms with both sig peak and trough 
        if pos_sigSTNpairs(i) == 1 & neg_sigSTNpairs(i) == 1
            pos_sigSTNpairs(i) = 0;
            neg_sigSTNpairs(i) = 0;
            STN_peak_and_trough = STN_peak_and_trough + 1;
        end
        % total CCFs
        total_STN_ccf = total_STN_ccf + STN_ccf{i};

    end
    
    for i = 1:n_GPe_pairs
        fprintf(1, 'GPe xcorr loop %d\n', i);
        
        [GPe_ccf{i},xGPe,f1,f2,n_pairs] = LIF_xcorr(GPe_times{GPe_pairs(i,1)},GPe_times{GPe_pairs(i,2)},binsize,[0 time_seconds],maxlag);
        total_GPe_pairs = total_GPe_pairs + n_pairs;
        [mGPe,sGPe,spGPe,stGPe] = xcorr_stats(GPe_ccf{i},n_pairs,p);
        
        % test for significance in cross-correlograms
        G_pos_sig_bins = find(GPe_ccf{i} >= spGPe);
        G_neg_sig_bins = find(GPe_ccf{i} < stGPe);
        
        if method1
            % METHOD 1: number of bins greater than confidence interval exceeds the expected amount
            if length(G_pos_sig_bins) > num_sig_bins 
                pos_sigGPepairs(i) = 1;   
            end
            if length(G_neg_sig_bins) > num_sig_bins 
                neg_sigGPepairs(i) = 1;   
            end
            
        else
            % METHOD 2: number of consecutive significant bins is greater than fixed threshold (3 - Raz et al 2000)
            for loop2 = 1:num_sig_bins
                G_pos_sig_bins = find(diff(G_pos_sig_bins) == 1); % indices of adjacent bins - reduce                
                G_neg_sig_bins = find(diff(G_neg_sig_bins) == 1); % indices of adjacent bins - reduce                
            end
            
            if ~isempty(G_pos_sig_bins)  pos_sigGPepairs(i) = 1; end 
            if ~isempty(G_neg_sig_bins)  neg_sigGPepairs(i) = 1; end
        end
       
        % omit correlograms with both sig peak and trough 
        if pos_sigGPepairs(i) == 1 & neg_sigGPepairs(i) == 1
            pos_sigGPepairs(i) = 0;
            neg_sigGPepairs(i) = 0;
            GPe_peak_and_trough = GPe_peak_and_trough + 1;
        end

        % total CCFs
        total_GPe_ccf = total_GPe_ccf + GPe_ccf{i};
    end
    
    
    mean_STN_ccf = total_STN_ccf ./ n_STN_pairs;
    mean_GPe_ccf = total_GPe_ccf ./ n_GPe_pairs;
    
    prop_POS_sig_STN_pairs = sum(pos_sigSTNpairs) / n_STN_pairs
    prop_NEG_sig_STN_pairs = sum(neg_sigSTNpairs) / n_STN_pairs
    prop_POS_sig_GPe_pairs = sum(pos_sigGPepairs) / n_GPe_pairs
    prop_NEG_sig_GPe_pairs = sum(neg_sigGPepairs) / n_GPe_pairs
    

    if do_total_xcorr_stats
        [mtotalSTN,stotalSTN,sptotalSTN,sttotalSTN] = xcorr_stats(total_STN_ccf,total_STN_pairs);
        [mtotalGPe,stotalGPe,sptotalGPe,sttotalGPe] = xcorr_stats(total_GPe_ccf,total_GPe_pairs);
        if do_display
            figure
            subplot(211),bar(xSTN,total_STN_ccf,1)
            axis([-0.1 0.1 min(total_STN_ccf)  max(total_STN_ccf)])
            line([-0.1 0.1],[sptotalSTN sptotalSTN],'Color',[1 0 0]);
            subplot(212),bar(xGPe,total_GPe_ccf,1)
            axis([-0.1 0.1 min(total_GPe_ccf)  max(total_GPe_ccf)])
            line([-0.1 0.1],[sptotalGPe sptotalGPe],'Color',[1 0 0]);
            drawnow
        end
    end
end

%%%%%%%%%%%%% do xcorr between nuclei %%%%%%%%%%%%%

pos_sigSTNGPepairs = zeros(n_inter_pairs,1);
neg_sigSTNGPepairs = zeros(n_inter_pairs,1);

if do_xcorr_inter
disp('doing xcorr - inter nucleus');
    STNGPe_peak_and_trough = 0; 

    for i = 1:n_inter_pairs
        [STNGPe_ccf{i},xSTNGPe,f1,f2,n_pairs] = LIF_xcorr(STN_times{STNGPe_pairs(i,1)},GPe_times{STNGPe_pairs(i,2)}, binsize, [0 time_seconds], maxlag);
        total_STNGPe_pairs = total_STNGPe_pairs + n_pairs;
        total_STNGPe_ccf = total_STNGPe_ccf + STNGPe_ccf{i};    
        [mSG,sSG,spSG,stSG] = xcorr_stats(STNGPe_ccf{i},n_pairs,p);

        % test for significance in cross-correlograms
        SG_pos_sig_bins = find(STNGPe_ccf{i} >= spSG);
        SG_neg_sig_bins = find(STNGPe_ccf{i} < stSG);
         if method1
            % METHOD 1: number of bins greater than confidence interval exceeds the expected amount
            if length(SG_pos_sig_bins) > num_sig_bins     % for just comparing all bins over correlogram
                pos_sigSTNGPepairs(i) = 1;    
            end
            if length(SG_neg_sig_bins) > num_sig_bins 
                neg_sigSTNGPepairs(i) = 1;   
            end            
        else
         % METHOD 2: number of consecutive significant bins is greater than fixed threshold (3 - Raz et al 2000)
            for loop2 = 1:num_sig_bins
                SG_pos_sig_bins = find(diff(SG_pos_sig_bins) == 1); % indices of adjacent bins - reduce                
                SG_neg_sig_bins = find(diff(SG_neg_sig_bins) == 1); % indices of adjacent bins - reduce      
                
            end
            
            if ~isempty(SG_pos_sig_bins)  pos_sigSTNGPepairs(i) = 1; end   % then there is at least one occurence of a num_sig_bins number of consecutive significant bins
            if ~isempty(SG_neg_sig_bins)  neg_sigSTNGPepairs(i) = 1; end
        end

        if pos_sigSTNGPepairs(i) == 1 & neg_sigSTNGPepairs(i) == 1
            pos_sigSTNGPepairs(i) = 0;
            neg_sigSTNGPepairs(i) = 0;
            STNGPe_peak_and_trough = STNGPe_peak_and_trough + 1;
        end

    end
    mean_STNGPe_ccf = total_STNGPe_ccf ./ n_inter_pairs;

    prop_POS_sig_STNGPe_pairs = sum(pos_sigSTNGPepairs) / n_inter_pairs
    prop_NEG_sig_STNGPe_pairs = sum(neg_sigSTNGPepairs) / n_inter_pairs


   
    if do_total_xcorr_stats
        [mSTNGPe,sSTNGPe,spSTNGPe,stSTNGPe] = xcorr_stats(total_STNGPe_ccf,total_STNGPe_pairs);
        if do_display
            figure
            bar(xSTNGPe,total_STNGPe_ccf,1)
            axis([-0.1 0.1 min(total_STNGPe_ccf)  max(total_STNGPe_ccf)])
            line([-0.1 0.1],[spSTNGPe spSTNGPe],'Color',[0 0 0]);
        end
    end
end

%%% do auto-correlograms %%%%
if do_acorr
    STN_acf = cell(n_STN,1);
    GPe_acf = cell(n_GPe,1);
    for i = 1:n_STN
        fprintf(1, 'STN autocorrelation loop count %d\n', i);
        [STN_acf{i},xSTNacf,f1,f2,n_pairs] = LIF_xcorr(STN_times{i},STN_times{i}, binsize, [0 time_seconds], maxlag);    
    end
    
    % [mSTN,sSTN,spSTN,stSTN] = xcorr_stats(STN_ccf{i},n_pairs,p);
    for i = 1:n_GPe
        fprintf(1, 'GPe autocorrelation loop count %d\n', i);
        [GPe_acf{i},xGPeacf,f1,f2,n_pairs] = LIF_xcorr(GPe_times{i},GPe_times{i}, binsize, [0 time_seconds], maxlag);    
    end
    
    % [mGPe,sGPe,spGPe,stGPe] = xcorr_stats(GPe_acf{i},n_pairs,p);    
end

%%%%%%%%%%%%%%% coherence analysis %%%%%%%%%%%%%
%%%% currently GPe and STN-GPe only %%%%%%%%%%%%%%%%%%%%%%
if do_coherence
   
    %%% do GPe coherence %%%
    GPe_cohere_res = struct('freqs', [], 'cohere', [], 'errs', [], 'peak_freq', [], 'peak_cohere', [], 'meanC', []);
        
    disp('doing GPe coherence');
    for i = 1:n_GPe_pairs
        fprintf(1, 'GPe coherence loop count %d\n', i);
        
        %% GPe %%% 
        data1 = GPe_data(GPe_pairs(i,1)).times;
        data2 = GPe_data(GPe_pairs(i,2)).times;
        [r1 c] = size(data1);
        [r2 c] = size(data2);
        r = min(r1,r2);
        if r >= Min_spikesC
            [C,phi,f,confC,phierr,Cerr]=coherencypt(data1, data2, tapersC, padC, FsC, freq_rangeC, err_barsC, 0,fscorrC);
            GPe_cohere_res(i).freqs = f;
            GPe_cohere_res(i).cohere = C;
            GPe_cohere_res(i).errs = Cerr;
            GPe_cohere_res(i).meanC = mean(C);

            GPe_err_lo = Cerr(1,:);
            sig_test_GPe_cohere = GPe_err_lo - mean(C);
            [GPe_peak GPe_peak_index] = max(sig_test_GPe_cohere);
            
            if GPe_peak >0
                GPe_cohere_res(i).peak_freq = f(GPe_peak_index);
                GPe_cohere_res(i).peak_cohere = GPe_peak;
            else
                GPe_cohere_res(i).peak_freq = NaN;
                GPe_cohere_res(i).peak_cohere = NaN;
            end
        else
            GPe_cohere_res(loop).freqs = [];
            GPe_cohere_res(loop).cohere = [];
            GPe_cohere_res(loop).errs = [];
            GPe_cohere_res(loop).meanC = [];
            GPe_cohere_res(loop).peak_freq = NaN;
            GPe_cohere_res(loop).peak_cohere = NaN;
        end
    end

    
    
    %%%% do mean GP Coherence %%%%%
    
    GPe_coheres = {GPe_cohere_res.cohere};
    GPe_Cfreqs = {GPe_cohere_res.freqs};
    N_freqs = cellfun('length', GPe_Cfreqs);
    N_non_null_pairsGP = length(N_freqs);
    [max_N_fs, index] = max(N_freqs);
    GPe_mu_cohere = zeros(max_N_fs,1);
    GPe_Cfreqs_base = GPe_Cfreqs{index};
    keyboard
    for loop=1:N_non_null_pairsGP
        if N_freqs(loop) > 0
            if N_freqs(loop) ~= max_N_fs
                Coheres = interp1(GPe_Cfreqs{loop}, GPe_coheres{loop}, GPe_Cfreqs_base);
                Coheres = Coheres';
            else
                Coheres = GPe_coheres{loop};
            end
        else
            Coheres = zeros(max_N_fs, 1);
        end
        GPe_mu_cohere = GPe_mu_cohere + Coheres;
    end
    GPe_mu_cohere = GPe_mu_cohere ./ N_non_null_pairsGP;
    
    sig_GPe_peak_Cfreqs = [GPe_cohere_res.peak_freq];
    GPe_sig_Cfs = sig_GPe_peak_Cfreqs(find(sig_GPe_peak_Cfreqs > 0));
    lowgc = GPe_sig_Cfs(GPe_sig_Cfs < 5);
    GPe_Cs = [GPe_cohere_res.meanC];

    %%%%% plot GPe coherence %%%%%
    if do_display
        
        figure
        subplot(2,1,1)
        hist(GPe_sig_Cfs, 20);
        title('significant maximum frequencies - coherence Gpe');
        xlabel('peak freq');

        subplot(2,1,2)
       
        hist(lowgc, 15);
        title('significant maximum frequencies in range 0-5Hz - coherence GPe');
        xlabel('peak freq');

        figure
	subplot(2,1,1)
        plot(GPe_Cfreqs_base, GPe_mu_cohere);
        title('Mean GPe coherence');
        xlabel('frequency');

        subplot(2,1,2)
        
        hist(GPe_Cs, 40);
	title('histogram of mean coherence values')
        drawnow
    end
    fprintf(1, 'mean GPe coherence %f\n\n', mean(GPe_Cs));

    %%%%% do GPe-STN coherence %%%%%

    STNGPe_cohere_res = struct('freqs', [], 'cohere', [], 'errs', [], 'peak_freq', [], 'peak_cohere', [], 'meanC', []);
   
    disp('doing STN-GPe coherence');
    for i = 1:n_inter_pairs
        fprintf(1, 'coherence loop count %d\n', i);
        
        %% STN-GPe %%% 
        data1 = GPe_data(STNGPe_pairs(i,2)).times;
        data2 = STN_data(STNGPe_pairs(i,1)).times;
        [r1 c] = size(data1);
        [r2 c] = size(data2);
        r = min(r1,r2);
        if r >= Min_spikesC
            [C,phi,f,confC,phierr,Cerr]=coherencypt(data1, data2, tapersC, padC, FsC, freq_rangeC, err_barsC, 0,fscorrC);
            STNGPe_cohere_res(i).freqs = f;
            STNGPe_cohere_res(i).cohere = C;
            STNGPe_cohere_res(i).errs = Cerr;
            STNGPe_cohere_res(i).meanC = mean(C);

            STNGPe_err_lo = Cerr(1,:);
            sig_test_STNGPe_cohere = STNGPe_err_lo - mean(C);
            [STNGPe_peak STNGPe_peak_index] = max(sig_test_STNGPe_cohere);
            
            if STNGPe_peak >0
                STNGPe_cohere_res(i).peak_freq = f(STNGPe_peak_index);
                STNGPe_cohere_res(i).peak_cohere = STNGPe_peak;
            else
                STNGPe_cohere_res(i).peak_freq = NaN;
                STNGPe_cohere_res(i).peak_cohere = NaN;
            end
        else
            STNGPe_cohere_res(loop).freqs = [];
            STNGPe_cohere_res(loop).cohere = [];
            STNGPe_cohere_res(loop).errs = [];
            STNGPe_cohere_res(loop).meanC = [];
            STNGPe_cohere_res(loop).peak_freq = NaN;
            STNGPe_cohere_res(loop).peak_cohere = NaN;
        end
    end

 %%%% do mean STN-GP Coherence %%%%%
    
    STNGPe_coheres = {STNGPe_cohere_res.cohere};
    STNGPe_Cfreqs = {STNGPe_cohere_res.freqs};
    N_freqs = cellfun('length', STNGPe_Cfreqs);
    N_non_null_pairsSTNGP = length(N_freqs);
    [max_N_fs, index] = max(N_freqs);
    STNGPe_mu_cohere = zeros(max_N_fs,1);
    STNGPe_Cfreqs_base = STNGPe_Cfreqs{index};
    for loop=1:N_non_null_pairsSTNGP
        if N_freqs(loop) > 0
            if N_freqs(loop) ~= max_N_fs
                Coheres = interp1(STNGPe_Cfreqs{loop}, STNGPe_coheres{loop}, STNGPe_Cfreqs_base);
                Coheres = Coheres';
            else
                Coheres = STNGPe_coheres{loop};
            end
        else
            Coheres = zeros(max_N_fs, 1);
        end
        STNGPe_mu_cohere = STNGPe_mu_cohere + Coheres;
    end
    STNGPe_mu_cohere = STNGPe_mu_cohere ./ N_non_null_pairsSTNGP;


%%%% plot STN GPe coherences %%%%%

    sig_STNGPe_peak_Cfreqs = [STNGPe_cohere_res.peak_freq];
    STNGPe_sig_Cfs = sig_STNGPe_peak_Cfreqs(find(sig_STNGPe_peak_Cfreqs > 0));
    lowsgc = STNGPe_sig_Cfs(STNGPe_sig_Cfs < 5);
    STN_GPe_Cs = [STNGPe_cohere_res.meanC];

    if do_display
        figure
        subplot(2,1,1)
        hist(STNGPe_sig_Cfs, 20);
        title('significant maximum  frequencies - coherence STN-Gpe');
        xlabel('peak freq');

        subplot(2,1,2)
        hist(lowsgc, 15);
        title('significant maximum  frequencies in range 0-5Hz - coherence STN-GPe');
        xlabel('peak freq');

        figure
	subplot(2,1,1)
        plot(STNGPe_Cfreqs_base, STNGPe_mu_cohere);
        title('Mean STN-GPe coherence')
        xlabel('frequency');

	subplot(2,1,2)
        
        hist(STN_GPe_Cs, 40);
        title('histogram of mean coherence values')
        drawnow
    end

    fprintf(1, 'mean STN-GPe coherence %f\n\n', mean(STN_GPe_Cs));

end

%%%%%%%%%%%%%%%%%%Spike triggered averaging %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if do_spike_trig_avg

    ctx_times = cell(neurons_per_nucleus,1);
    ctx_IFRbins = cell(neurons_per_nucleus,1);
    ctx_mean = [];
    
    % generate pseudo-EEG using mean of smoothed single-neuron
    % firing rate estimates
    fprintf(1, 'Computing pseudo cortical EEG \n');
    for loop = 1:neurons_per_nucleus
        ctx_times{loop} = t_in_t(in_n == EXT(loop));
        [ctx_IFRbins{loop},ctx_timepoints] = LIF_firingrate(ctx_times{loop}',smooth_win_size,time_seconds,dt,step_size,window);
        if loop==1 ctx_mean = ctx_IFRbins{loop} ./ neurons_per_nucleus; 
        else
            ctx_mean = ctx_mean + ctx_IFRbins{loop} ./ neurons_per_nucleus;  
        end
    end
   

    smooth_steps = length(ctx_timepoints);
    new_dt = step_size * dt;
    win_in_steps = trig_win_half ./ new_dt;
    
    avg_x_axis = linspace(-win_in_steps,win_in_steps,win_in_steps*2+1) .* new_dt;
    
    GPe_avg_wvfrm = zeros(n_GPe,win_in_steps*2+1);
    STN_avg_wvfrm = zeros(n_STN,win_in_steps*2+1);
     

    for loop1 = 1:n_GPe
        GPe_t_stamps = GPe_times{loop1};
        GPe_num_spikes = length(GPe_t_stamps);
          
        
        for loop2 = 1:GPe_num_spikes
            start_t = round(GPe_t_stamps(loop2) ./ new_dt) - win_in_steps;
            end_t = round(GPe_t_stamps(loop2) ./ new_dt) + win_in_steps;
            if start_t > 0 & end_t <= smooth_steps 
                GPe_avg_wvfrm(loop1,:) = GPe_avg_wvfrm(loop1,:) + ctx_mean(start_t:end_t);
            end
        end
        GPe_avg_wvfrm(loop1,:) = GPe_avg_wvfrm(loop1,:) ./ GPe_num_spikes;
    end
    
    for loop1 = 1:n_STN   
        STN_t_stamps = STN_times{loop1};
        STN_num_spikes = length(STN_t_stamps);

        for loop2 = 1:STN_num_spikes
            start_t = round(STN_t_stamps(loop2) ./ new_dt) - win_in_steps;
            end_t = round(STN_t_stamps(loop2) ./ new_dt) + win_in_steps;
            if start_t > 0 & end_t <= smooth_steps 
                STN_avg_wvfrm(loop1,:) = STN_avg_wvfrm(loop1,:) + ctx_mean(start_t:end_t);
            end
        end
        STN_avg_wvfrm(loop1,:) = STN_avg_wvfrm(loop1,:) ./ STN_num_spikes;
    end

%     figure
%     plot(avg_x_axis,GPe_avg_wvfrm(2,:))
%     axis([-win_in_steps*new_dt win_in_steps*new_dt 0 max(ctx_mean)]);
%     drawnow
%     disp('press to continue')
%     pause;
end


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

analysis_modes = struct('means', do_means, 'isi', do_isi, 'spectra', do_spectra,...
        'bursts', do_bursts, 'xcorr_intra', do_xcorr_intra, 'xcorr_inter', do_xcorr_inter,...
        'coherence',do_coherence,'spike_trig_avg',do_spike_trig_avg);

an_fname = [pathroot an_fname];

if save_full_analysis
    save(an_fname, 'analysis_modes','t_in_t','ISIbins','limits',...
        'freq_range', 'tapers', 'pad', 'Fs', 'err_bars', 'fscorr', 'binsize','maxlag','p','n_STN','n_GPe','n_GPi');
    if do_means
        save(an_fname, 'mean_STN','mean_GPe','mean_GPi','STN_Hz','GPe_Hz','GPi_Hz',...
            'sem_STN', 'std_STN', 'sem_GPe', 'std_GPe', 'sem_GPi', 'std_GPi', 'CV_STN', 'CV_GPe', 'CV_GPi', '-append');
    end
    if do_isi
        save(an_fname, 'STN_rates','STN_ISIhist','STN_x_isi','STN_x_hist',...
            'GPe_rates','GPe_ISIhist','GPe_x_isi','GPe_x_hist', '-append');
    end
    if do_spectra
        save(an_fname, 'STN_spect_res', 'GPe_spect_res',...
             'STN_sig_fs', 'GPe_sig_fs','-append');
    end
    if do_bursts
        save(an_fname, 'N_STN_bursts', 'N_GPe_bursts',...
            'STN_bursts','STN_burst_isis', 'GPe_bursts','GPe_burst_isis', '-append');
    end
    if do_xcorr_intra
        save(an_fname, 'n_intra_pairs', 'STN_pairs','GPe_pairs', ...
            'STN_ccf','xSTN','mean_STN_ccf','GPe_ccf','xGPe','mean_GPe_ccf',...
            'prop_POS_sig_STN_pairs','prop_NEG_sig_STN_pairs','prop_POS_sig_GPe_pairs','prop_NEG_sig_GPe_pairs',...
            'STN_peak_and_trough','GPe_peak_and_trough', '-append');
    end
    if do_xcorr_inter
        save(an_fname,'n_inter_pairs','STNGPe_pairs',...
            'STNGPe_ccf','xSTNGPe','mean_STNGPe_ccf',... 
            'prop_POS_sig_STNGPe_pairs','prop_NEG_sig_STNGPe_pairs','STNGPe_peak_and_trough','-append');
    end
    if do_acorr 
        save(an_fname,'STN_acf','xSTNacf','GPe_acf','xGPeacf','-append');
    end
    if do_coherence
        save(an_fname, 'n_intra_pairs','GPe_pairs', 'GPe_cohere_res', 'GPe_Cfreqs_base', 'GPe_mu_cohere',...
	     'n_inter_pairs','STNGPe_pairs', 'STNGPe_cohere_res', 'STNGPe_Cfreqs_base', 'STNGPe_mu_cohere', '-append');
    end
   if do_spike_trig_avg
        save(an_fname, 'ctx_mean', 'GPe_avg_wvfrm','STN_avg_wvfrm', '-append');
    end
end
if save_summary_analysis
    
    save(an_fname, 'neurons_per_nucleus', 'maxlag','analysis_modes');

    if do_means
        save(an_fname, 'mean_STN','mean_GPe','mean_GPi','STN_Hz', 'sem_STN', 'std_STN',...
            'GPe_Hz', 'sem_GPe', 'std_GPe', 'GPi_Hz', 'sem_GPi', 'std_GPi', 'CV_STN', 'CV_GPe', 'CV_GPi', '-append');
    end
    if do_isi
        save(an_fname, 'STN_times', 'GPe_times', 'rs_STN', 'rs_GPe', 'rs_in', 'STN_rates', 'STN_x_isi',...
	     'GPe_rates', 'GPe_x_isi', '-append');
    end
    if do_spectra
        save(an_fname, 'STN_freqs_base', 'STN_mu_power', 'GPe_freqs_base', 'GPe_mu_power',...
            'STN_sig_fs', 'GPe_sig_fs', 'lows', 'lowg', '-append');
    end
    if do_bursts
        save(an_fname, 'N_STN_bursts', 'N_GPe_bursts', '-append');
    end
    if do_xcorr_intra 
        save(an_fname,'xSTN', 'total_STN_ccf', 'xGPe', 'total_GPe_ccf', 'n_intra_pairs',...
            'prop_POS_sig_STN_pairs','prop_NEG_sig_STN_pairs','prop_POS_sig_GPe_pairs','prop_NEG_sig_GPe_pairs',...
            'STN_peak_and_trough','GPe_peak_and_trough','-append');
    end
    if do_xcorr_inter
        save(an_fname,'n_inter_pairs', 'total_STNGPe_ccf','xSTNGPe',... 
            'prop_POS_sig_STNGPe_pairs','prop_NEG_sig_STNGPe_pairs','STNGPe_peak_and_trough', '-append');
    end
    if do_acorr 
        save(an_fname,'STN_acf','xSTNacf','GPe_acf','xGPeacf','-append');
    end
    if do_coherence
        save(an_fname, 'n_intra_pairs', 'N_non_null_pairsGP','N_non_null_pairsSTNGP',  'GPe_Cfreqs_base', 'GPe_mu_cohere', 'GPe_Cs', 'STN_GPe_Cs',...
	      'STNGPe_Cfreqs_base', 'STNGPe_mu_cohere', 'GPe_sig_Cfs', 'STNGPe_sig_Cfs',...
	      'lowgc', 'lowsgc', '-append');
    end
    if do_spike_trig_avg
        save(an_fname, 'ctx_mean', 'GPe_avg_wvfrm', 'STN_avg_wvfrm', '-append');
    end
   
end