clear all;
close all;

fname = '2013_02_26_Yam-22565_no-drug_01_nowaveform';
fid_label = fname;
dirname = 'C:\Users\USER\Documents\mingfaifong\opticalstim ctrl cnqx apv\matfiles\';
% timestamps = [3*3600 15*3600];
load([dirname fname]);
timestamps = [0 max(time)]


%%
clearvars -except type dirID fidlist time timestamps channel excludechannels thresASDR fid_label unit meta
savedir = 'C:\Users\USER\Documents\mingfaifong\opticalstim ctrl cnqx apv\burst analysis\';
excludechannels = [1 8 33 57 64];
channel = channel(ismember(channel,excludechannels)==0);
time = time(ismember(channel,excludechannels)==0);

%% 
dt = 0.01; %1s

dat.time = time(time>=timestamps(1) & time<timestamps(2));
dat.channel = channel(time>=timestamps(1) & time<timestamps(2));
dat.unit= unit(time>=timestamps(1) & time<timestamps(2));
dat.bin = timestamps(1):dt:timestamps(2);
dat.asdr = histc(dat.time,dat.bin)/dt;

type = input('treatment: ');

datstarttime = timestamps(1);
datendtime = timestamps(2);

close all;
figure(1);
stairs(dat.bin+dt/2,dat.asdr,'-k');
xlim([0 datendtime])

thresASDR = input('ASDR threshold: ');

%%
% thresASDR = 1000;
minIBI = 2;

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

hold on;
plot(dat.burst.time,dat.asdr(dat.burst.ind),'*b');


set(gcf,'Position',[4 549 2547 373])
xlim([0 datendtime])
set(gca,'FontSize',12,'XTick',[0:1:24]*3600,'XTickLabel',[0:1:24],'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');

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

%% get burst snippet times

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

figure(5); 
dat.burst.snippet.FR = [];
dat.burst.snippet.spiketimes = [];
dat.burst.snippet.sniptime = 0;
subplot(1,2,1); hold off;
subplot(1,2,2); hold off;
burstcolor = jet(length(dat.burst.ind));

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

set(gca,'FontSize',12,'TickDir','out','TickLength',[0.0025 0.025],'XTick',[snip_min:1:snip_max],'yscale','log');
box off;
grid on;
xlim([snip_min snip_max]);

subplot(1,2,1);
set(gca,'FontSize',12,'TickDir','out','TickLength',[0.0025 0.025],'XTick',[snip_min:1:snip_max]);
box off;
grid on;
xlim([snip_min snip_max]);
xlabel('Time (sec)','FontSize',16);
ylabel('ASDR (Hz)','FontSize',16);

title([fid_label ' ' type ', Burst Overlay'],'Interpreter','none');

set(gcf,'Position',[1281 300 1280 957/2]);

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

%% interburst interval histogram
figure;

dat.IBI = diff(dat.burst.time);
xmax = ceil(max(dat.IBI)/100)*100;
xint = ceil(xmax/10/60); %minutes

dtIBI = 20;
binIBI = 0:dtIBI:max(dat.IBI);

hist(dat.IBI,binIBI);
h = findobj(gca,'Type','patch');
set(h,'FaceColor','b','EdgeColor','k','LineWidth',1)
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;
xlabel('Interburst Interval (min)','FontSize',16);
ylabel('Counts','FontSize',16);
title([fid_label ' IBI histogram, with dat (24 hours)'],'Interpreter','none');
% xlim([0 10*60])

set(gcf,'Position',[360 388 589 534]);

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

%% interburst interval over time
figure(3);
yint = xint;
ymax = xmax;

hold off;
plot(dat.burst.time(2:end),dat.IBI,'bo');
xlim([0 datendtime])
ylim([0 ymax]);
set(gca,'FontSize',12,'XTick',[0 datendtime],'XTickLabel',[0 datendtime]/3600, ...
    'YTick',[0:60*yint:ymax],'YTickLabel',[0:yint:ymax/60], ...
    'TickDir','out','TickLength',[0.0025 0.025]);
box off;
grid on;
axis square;
xlabel('Time (hours)','FontSize',16);
ylabel('Interburst Interval (min)','FontSize',16);
set(gcf,'Position',[360 388 589 534]);
% ylim([0 70*60]);

title([fid_label ', IBIs over time'],'Interpreter','none');

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


% %% interburst interval FR
% for i = 1:size(dat.burst.snippet.timestamps,1)
%     if i == 1
%         dat.IBI.bounds(i,1) = datstarttime;
%     else
%         dat.IBI.bounds(i,1) = dat.burst.snippet.timestamps(i,2);
%     end
%     if i == size(dat.burst.snippet.timestamps,1)
%         dat.IBI.bounds(i,2) = datendtime;        
%     else
%         dat.IBI.bounds(i,2) = dat.burst.snippet.timestamps(i+1,1);
%     end
% end
% 
% % if strcmp(type,'CNQXstim1X')
% %     bin_sec = [-6*3600:3*3600:0,300:3*3600:24*3600+300];
% % else
%     bin_sec = [0:3*3600:datendtime];
% %     bin_sec = [-6*3600:3*3600:24*3600];
% % end
% 
% IBI.bounds = [pre.IBI.bounds; dat.IBI.bounds];
% 
% for i = 1:length(bin_sec)
%     for j = 1:size(IBI.bounds,1)
%         if bin_sec(i)>IBI.bounds(j,1) && bin_sec(i)<IBI.bounds(j,2)
%             IBI.bounds(j+2:end+1,:) = IBI.bounds(j+1:end,:);
%             IBI.bounds(j:j+1,:) = [IBI.bounds(j,1),bin_sec(i);bin_sec(i),IBI.bounds(j,2)];
%         end
%     end
% end
% 
% m = 1;
% dur_IBI = zeros(1,length(bin_sec)-1);
% count_spkIBI = zeros(1,length(bin_sec)-1);
% for n = 1:size(IBI.bounds,1)
%     if diff(IBI.bounds(n,:))<0
%         continue;
%     end
%     if IBI.bounds(n,1)<bin_sec(m)
%         continue;
%     elseif IBI.bounds(n,2)>bin_sec(m+1)
%         m = m+1;
%     end
%     dur_IBI(1,m) = dur_IBI(1,m) + diff(IBI.bounds(n,:));
%     count_spkIBI(1,m) = count_spkIBI(1,m) + length(time(time>=IBI.bounds(n,1) & time<IBI.bounds(n,2)));
% end
% 
% if strcmp(type,'CNQXstim1X')
%     dur_IBI = [dur_IBI(1:2) dur_IBI(4:end)];
%     count_spkIBI = [count_spkIBI(1:2) count_spkIBI(4:end)];
% end
% 
% 
% %% spk an burst counts
% count_spk = [length(pre.time(pre.time>=-6*3600 & pre.time<-3*3600)),...
%     length(pre.time(pre.time>=-3*3600 & pre.time<0)),...
%     length(dat.time(dat.time>=datstarttime & dat.time<datstarttime+3*3600)),...
%     length(dat.time(dat.time>=datstarttime+3*3600 & dat.time<datstarttime+6*3600)),...
%     length(dat.time(dat.time>=datstarttime+6*3600 & dat.time<datstarttime+9*3600)),...
%     length(dat.time(dat.time>=datstarttime+9*3600 & dat.time<datstarttime+12*3600)),...
%     length(dat.time(dat.time>=datstarttime+12*3600 & dat.time<datstarttime+15*3600)),...
%     length(dat.time(dat.time>=datstarttime+15*3600 & dat.time<datstarttime+18*3600)),...
%     length(dat.time(dat.time>=datstarttime+18*3600 & dat.time<datstarttime+21*3600)),...
%     length(dat.time(dat.time>=datstarttime+21*3600 & dat.time<datstarttime+24*3600))];
% 
% count_burst = [length(pre.burst.time(pre.burst.time>=-6*3600 & pre.burst.time<-3*3600)),...
%     length(pre.burst.time(pre.burst.time>=-3*3600 & pre.burst.time<0)),...
%     length(dat.burst.time(dat.burst.time>=datstarttime & dat.burst.time<datstarttime+3*3600)),...
%     length(dat.burst.time(dat.burst.time>=datstarttime+3*3600 & dat.burst.time<datstarttime+6*3600)),...
%     length(dat.burst.time(dat.burst.time>=datstarttime+6*3600 & dat.burst.time<datstarttime+9*3600)),...
%     length(dat.burst.time(dat.burst.time>=datstarttime+9*3600 & dat.burst.time<datstarttime+12*3600)),...
%     length(dat.burst.time(dat.burst.time>=datstarttime+12*3600 & dat.burst.time<datstarttime+15*3600)),...
%     length(dat.burst.time(dat.burst.time>=datstarttime+15*3600 & dat.burst.time<datstarttime+18*3600)),...
%     length(dat.burst.time(dat.burst.time>=datstarttime+18*3600 & dat.burst.time<datstarttime+21*3600)),...
%     length(dat.burst.time(dat.burst.time>=datstarttime+21*3600 & dat.burst.time<datstarttime+24*3600))];
% z = [count_spk,count_burst,count_spkIBI,dur_IBI];

%%
% close all;
save([savedir fid_label '_burst.mat']);