clear all;
close all;

fid = 'quimby22930_15-19DIV';
i_well = 6; %A=1, B=2, C=3, D=4, E=5, F=6
fid_label = [fid '_w' num2str(i_well)];
load([fid '.mat']);

w(1,:) = [3 4 5 12 13 20 21 28 29];
w(2,:) = [7 14 15 16 22 23 24 30 31];
w(3,:) = [38 39 40 46 47 48 55 56 63];
w(4,:) = [36 37 44 45 52 53 60 61 62];
w(5,:) = [34 35 41 42 43 49 50 51 58];
w(6,:) = [2 9 10 17 18 19 25 26 27];

time = time(ismember(channel,w(i_well,:))==1);
channel = channel(ismember(channel,w(i_well,:))==1);

xlim_plots = [-12*3600 36*3600];

%% Determining Burst Threshold
dt = 0.01; % 10ms bin
bin = floor(min(time)/dt)*dt:dt:ceil(max(time)/dt)*dt;
asdr = histc(time,bin)/dt;

close all;
figure(1);
set(gcf,'Position',[4 549 2547 373])
stairs(bin+dt/2,asdr,'-k');
xlim(xlim_plots);

thresASDR = input('ASDR threshold: ');
minIBI = 30; %seconds

%% Finding/Plotting Burst Times

x = find(asdr>thresASDR);
burst.ind = [];
burst.ind = x(1);
c = 1;
for i = 1:length(x)
    if bin(x(i))-minIBI>bin(burst.ind(c))
        c = c+1;
        burst.ind(c) = x(i);
    end
end;
burst.time = bin(burst.ind);

% figure(1);
% hold on;
% plot(burst.time,asdr(burst.ind),'*r');
% set(gcf,'Position',[4 549 2547 373])
% xlim(xlim_plots);
% set(gca,'FontSize',12,...
%     'XTick',xlim_plots(1):3*3600:xlim_plots(2),...
%     'XTickLabel',[xlim_plots(1):3*3600:xlim_plots(2)]/3600,...
%     'TickDir','out','TickLength',[0.0025 0.025]);
% box off;
% xlabel('Time (hours)','FontSize',16);
% ylabel('Array-Wide Firing Rate (Hz)','FontSize',16);
% 
% title([fid_label ' ASDR (bin size=' num2str(dt) 's) w/ detected bursts (thres=' num2str(thresASDR) 'Hz)'],'Interpreter','none');

%% Get Burst Snippet Times

snip_min = -2;
snip_max = 8;
snippetWidth = [snip_min:dt:snip_max];

burst.snippet.FR = [];
burst.snippet.spiketimes = [];
burst.snippet.sniptime = 0;

for i = 1:length(burst.ind)
    if burst.ind(i)+snip_min/dt < 1
        burst.snippet.FR(i,:) = nan(size(snippetWidth));
        burst.snippet.timestamps(i,:) = [bin(1),bin(burst.ind(i)+snip_max/dt)];
    elseif burst.ind(i)+snip_max/dt > length(asdr)
        burst.snippet.FR(i,:) = nan(size(snippetWidth));
        burst.snippet.timestamps(i,:) = [bin(burst.ind(i)+snip_min/dt),bin(end)];
    else
        burst.snippet.FR(i,:) = asdr(burst.ind(i)+snip_min/dt:burst.ind(i)+snip_max/dt);
        burst.snippet.timestamps(i,:) = [bin(burst.ind(i)+snip_min/dt),bin(burst.ind(i)+snip_max/dt)];
    end
    burst.snippet.spiketimes = vertcat(burst.snippet.spiketimes,...
        time(time>burst.snippet.timestamps(i,1)...
        & time<burst.snippet.timestamps(i,2)));
    if i>1 && burst.snippet.timestamps(i,1)<burst.snippet.timestamps(i-1,2)
        extratime = burst.snippet.timestamps(i,2)-burst.snippet.timestamps(i-1,2);
    else
        extratime = 10;
    end
    burst.snippet.sniptime = burst.snippet.sniptime + extratime;
end

% figure(2);
% subplot(3,1,1);
% plot(snippetWidth,burst.snippet.FR(burst.time<0,:));
% subplot(3,1,2);
% plot(snippetWidth,burst.snippet.FR(burst.time>0 & burst.time<24*3600,:));
% ylabel('Array-wide Firing Rate (Hz)','FontSize',16)
% subplot(3,1,3);
% if(isempty(burst.snippet.FR(burst.time>24*3600,:))==0)
%     plot(snippetWidth,burst.snippet.FR(burst.time>24*3600,:));
% end
% xlabel('Time (sec)','FontSize',16);
% 
% for i = 1:3
%     subplot(3,1,i)
%     xlim([-0.5 1.5]);
%     set(gca,'FontSize',12,'TickDir','out','TickLength',[0.0025 0.025]);
%     box off;
% end

%% Interburst Interval Histogram
pre.IBI = diff(burst.time(burst.time<0));
treat.IBI = diff(burst.time(burst.time>=0 & burst.time<24*3600));
post.IBI = diff(burst.time(burst.time>=24*3600));
xmax = ceil(max([pre.IBI,treat.IBI,post.IBI])/100)*100;
xmax = 20*60;
xint = ceil(xmax/10/60); %minutes

dtIBI = 5; %minutes
binIBI = 0:dtIBI:max([pre.IBI,treat.IBI,post.IBI]);

% figure(3);
% 
% subplot(3,1,1);
% hist(pre.IBI,binIBI);
% h = findobj(gca,'Type','patch');
% set(h,'FaceColor','k','EdgeColor','k','LineWidth',1)
% title([fid_label ' IBI histogram, pre-treatment'],'Interpreter','none');
% 
% subplot(3,1,2);
% hist(treat.IBI,binIBI);
% h = findobj(gca,'Type','patch');
% set(h,'FaceColor','r','EdgeColor','k','LineWidth',1)
% title([fid_label ' IBI histogram, during treatment'],'Interpreter','none');
% 
% subplot(3,1,3);
% hist(post.IBI,binIBI);
% h = findobj(gca,'Type','patch');
% set(h,'FaceColor','b','EdgeColor','k','LineWidth',1)
% title([fid_label ' IBI histogram, post-treatment'],'Interpreter','none');
% xlabel('Interburst Interval (min)','FontSize',16);
%     
% for i = 1:3
%     subplot(3,1,i)
%     xlim([0 xmax]);
%     set(gca,'FontSize',12,...
%         'XTick',[0:60*xint:xmax],...
%         'XTickLabel',[0:1*xint:xmax/60],...
%         'TickDir','out','TickLength',[0.0025 0.025]);
%     box off;
%     ylabel('Counts','FontSize',16);
% end
% 
% set(gcf,'Position',[360 388 589 534]);

%% Counting Spikes and Bursts
inc = 1*3600; %in seconds
bin_count = floor(min(time)/inc)*inc:inc:ceil(max(time)/inc)*inc;

spk_count = histc(time,bin_count);
burst_count = histc(burst.time,bin_count);

figure(4);
subplot(2,1,1);
stairs(bin_count/3600,spk_count,'LineWidth',2);
ylabel('# Spikes','FontSize',16);
subplot(2,1,2);
stairs(bin_count/3600,burst_count,'LineWidth',2)
ylabel('# Bursts','FontSize',16);
xlabel('Time (hours)','FontSize',16);

for i = 1:2
    subplot(2,1,i)
    set(gca,'FontSize',12,'TickDir','out','TickLength',[0.0025 0.025]);
    box off;
    grid on;
    xlim([bin_count(1),bin_count(end)]/3600);
end

spk_count = spk_count';
set(gcf,'Position',[13 307 1257 615]);

%% Save Burst Info
savedir = 'C:\Users\USER\Documents\CandaceLaw\analysis_ming\matfiles_findBursts\';
save([savedir fid_label '_burst'],'asdr','bin','binIBI','burst','channel','dt','dtIBI','extratime','fid','fid_label','i_well','minIBI',...
    'snip_max','snip_min','snippetWidth','thresASDR','time','w','xint','xlim_plots','xmax','inc','bin_count','spk_count','burst_count');

%% Save Figures

% savedir = 'C:\Users\USER\Documents\CandaceLaw\analysis_ming\figures_findBursts\';
% figure(1);
% export_fig([savedir fid_label '_ASDRwithBursts'],'-png','-transparent');
% saveas(gca,[savedir fid_label '_ASDRwithBursts'],'fig');
% 
% figure(2);
% export_fig([savedir fid_label '_BurstOverlay'],'-png','-transparent');
% saveas(gca,[savedir fid_label '_BurstOverlay'],'fig');
% 
% figure(3);
% export_fig([savedir fid_label '_HistogramIBI'],'-png','-transparent');
% saveas(gca,[savedir fid_label '_HistogramIBI'],'fig');

figure(4);
export_fig([savedir fid_label '_SpkBurstCountsOverTime'],'-png','-transparent');
saveas(gca,[savedir fid_label '_SpkBurstCountsOverTime'],'fig');