function [boot_stats, boot_results, correlations] = bootstrap(data,varargin)
rng('shuffle');
%%
% default values
do_plot=0;
bin_size=60; %sec
numOfIterations = 10000;
time_seg = 1200;
bsl = 0;
hour = data.networkFullHours;
%%
pvpmod(varargin);
hour = hour*1200;

dist = @(x,y) x-y; 
bsl_start = (bsl)*time_seg;
bsl_end = (bsl+1)*time_seg;
last_start = hour-time_seg;
last_end = hour;


for u=1:numel(data.unitSpikeTimes)
[r,c]=find(bsl_start<data.unitSpikeTimes{u} & data.unitSpikeTimes{u}<bsl_end);
bsl_spikes=data.unitSpikeTimes{u}(r)-bsl_start;
[r,c]=find(last_start<data.unitSpikeTimes{u} & data.unitSpikeTimes{u}<last_end);
last_spikes=data.unitSpikeTimes{u}(r)-last_start;

%%
bins=0:bin_size:time_seg;
bsl_freq=histc(bsl_spikes,bins)/bin_size;
bsl_freq=bsl_freq(1:numel(bins)-1);
end_freq=histc(last_spikes,bins)/bin_size;
end_freq=end_freq(1:numel(bins)-1);
if size(end_freq,2)>2
    end_freq=end_freq';
    end_freq=end_freq(1:numel(bins)-1);
end
if size(bsl_freq,2)>2
bsl_freq=bsl_freq';
bsl_freq=bsl_freq(1:numel(bins)-1);
end
mean_freq(u,1:4)=[mean(bsl_freq), std(bsl_freq), mean(end_freq), std(end_freq)];
real_dist= dist (mean(bsl_freq), mean(end_freq));
combined_freq= [bsl_freq, end_freq];
all_real_dist(u,1)=real_dist;
all_bsl_freq (u,1) = mean(bsl_freq);
all_end_freq (u,1) = mean(end_freq);
all_bsl_freq (u,2) = std(bsl_freq);
all_end_freq (u,2) = std(end_freq);
all_dist = zeros(1, numOfIterations);

%%
for i=1:numOfIterations
    tmp = randperm(size(bsl_freq,1)*2);
    firstPositions = tmp(1:size(bsl_freq,1));
    secondPositions = tmp(size(bsl_freq,1)+1:end);
    
    nonreal_bsl = combined_freq(firstPositions);
    nonreal_end = combined_freq(secondPositions);
    
    current_dist = dist(mean(nonreal_bsl) , mean(nonreal_end));
    all_dist(i) = current_dist;
    
end

confidence_interval = icdf('Normal',[0.025 0.975], mean(all_dist),std(all_dist) );

if do_plot==1
    figure
    hist(all_dist,100)
    line ([confidence_interval(1) confidence_interval(1)], [0 300], 'color', 'r')
    line ([confidence_interval(2) confidence_interval(2)], [0 300], 'color', 'r')
    line ([real_dist real_dist], [0 300], 'color', 'g')
end

if real_dist > confidence_interval(2)
    results(u,1) = 1; % significantly lower
elseif real_dist < confidence_interval(1)
    results(u,1) = -1; % significantly higher
else
    results(u,1) = 0; %not changed
end

%%




end

boot_results=[all_bsl_freq, all_end_freq, all_real_dist, results];
increased =  numel(find(results == -1))/numel(results);
decreased =  numel(find(results == 1))/numel(results);
no_change =  numel(find(results == 0))/numel(results);
boot_stats = [increased, no_change, decreased];
percent_change = 100* all_end_freq(:,1) ./ all_bsl_freq(:,1);
correlations= [all_bsl_freq(:,1) all_end_freq(:,1) percent_change(:,1) data.allUnitPC(:,3) data.allunitPeakValley(:,2) data.allunitPeakValley(:,3)];
end
